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FOREWORD 


This report, prepared by the Dynamics and Loads Section, Martin 
Marietta Corporation, Denver Division, under Contract NAS5-11996, 
presents the results of a study whose purpose was to develop a 
computer program system for dynamic simulation and stability 
analysis of passive and actively controlled spacecraft. The 
study was performed from May 1973 to April 1975 and was admin- 
istered by the National Aeronautics and Space Administration, 
Goddard Space Flight Center, Greenbelt, Maryland, under the di- 
rection of Mr. Joseph P, Young. 

The report is published in four volumes: 

Volume I - Theory 
Volume II - Program Users ^ Guide 
Volume III - Demonstration Problems 
Volume IV - Program Listing 
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ABSTRACT 


A theoretical development and associated digital computer program 
system for the dynamic simulation and stability analysis of pas- 
sive and actively controlled spacecraft is presented. The dynamic 
system (spacecraft) is modeled as an assembly of rigid and/or 
flexible bodies not necessarily in a topological tree configura- 
tion. The computer program system may be used to investigate 
total system dynamic characteristics including interaction effects 
between rigid and/or flexible bodies, control systems, and a 
wide range of environmental loadings. Additionally, the program 
system may be used for design of attitude control systems and for 
evaluation of total dynamic system performance including time do- 
main response and frequency domain stability analyses. 

Volume I presents the theoretical developments Including a des- 
cription of the physical system, the equations of dynamic equi- 
librium, discussion of kinematics and system topology, a complete 
treatment of momentum wheel coupling, and a discussion of gravity 
gradient and environmental effects. 

The development of synthesis and analysis techniques for the 
linearized system includes a discussion of the numerical linear- 
ization technique, procedures for definition of system transfer 
functions, and linear time domain response. 

Volume II is a program users* guide and includes a description of 
the overall digital program code, individual subroutines and a 
description of required program input and generated program out- 
put. 

Volume III presents the results of selected demonstration prob- 
lems that illustrate all program system capabilities. 

Volume IV contains a listing of the digital code. 
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Applications of such methods and program systems are numerous 
and include simulation of the Space Shuttle payload deplo3rment/ 
retrieval mechanism, solar panel array deployment, antenna de- 
ployment, analysis of raultlspin satellites, and analysis of 
large, highly flexible satellites. 

Our approach provides a general-purpose modeling capability for 
dynamic simulation and stability analysis of passive and actively 
controlled spacecraft. In particular, the following items are 
considered: (1) time domain solution of the nonlinear differen- 

tial equations of motion that describe total or nominal response* 
of the complete spacecraft system idealized as a collection of 
interconnected flexible (or rigid) bodies, (2) linearization of 
the governing equations by numerical means, (3) time domain so- 
lution of the linearized equations that describe the perturbation 
response of the complete spacecraft system about some predeter- 
mined (calculated or user-specified) nominal motion, C4) general 
frequency domain stability analysis corresponding to the linear- 
ized spacecraft representation, and (5) provision for arbitrary 
(explicitly time-dependent) loadings and environment interaction 
such as gravity gradient and thermally induced deformations re- 
sulting from solar radiation. 


B. DESCRIPTION OF THE PHYSICAL SYSTEM 


The physical system undergoing analysis may be generally descri- 
bed as a cluster of co ntiguous > flexible structures ^ (bodies) that 
comprise a mechanical system such as a spacecraft. The entire 
system (spacecraft) or portions thereof may be spinning or non- 
spinning. Member bodies of the spacecraft are capable of under- 
going large relative excursions such as those of appendage de- 
ployment, or rotor/stator motions. The general system of bodies 
is, by its inherent nature, a feedback system wherein inertial 
forces (such as those due to centrifugal and Coriolis accelera- 
tion) and the restoring and damping forces are motion dependent. 
Also, the system may possess a control system, wherein certain 
position and rate errors are actively controlled through use of 
reaction control jets, servomotors, or momentum wheels. 


* The total response of the dynamic system may be, in certain 
cases, considered to be equilibrium state motion (nominal 
response) plus perturbation motion with respect to the 
equilibrium state. 
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Bodies of the system may be interconnected by linear or nonlinear 
springs and dampers; they may be interconnected via a mechanism 
that is a combination of gimbal and slider block, or any combi 
nation of the above. Also, any two bodies of the system may be 
free (one from the other) and possess six degrees of relative 
motion freedom. Also, several or all of the six. degrees of rel- 
ative motion freedom, between two bodies, may be a prescribed 
function of time (including zero motion) . 

For purposes of further Introduction of the physical system, let 
us consider an illustrative example, such as the system of bodies 
of Figure I-l, and let this example be the prototype for all sub- 
sequent discussion and development . 

In Figure I-l, we have deliberately indicated a nontopologlcal 
tree configuration. There are five hinges and four bodies, thus 
one closed path. Consecutive integer labels are used for body 
reference points, for hinges, for sensor points, and for momen- 
tum wheels. The numerical orde r within each of the four label 
sets may be rando i^ However” it is understood that body _1 (which 
may be any body of the system) i s ass ocia ted with hinge 1 . 

For each body of the system, there is a body— fixed, right-handed 
reference frame, whose origin may be at the body's mass center 
or at some structural hard point on the body (a body's elastic 
deformation is measured in its reference frame). 


In this work a hinge is defined to be a pair of structural hard 
points (see Figure 1—2) with a point situated on each of two con- 
tiguous bodies. In Figure 1-2, point p and point q comprise a 
hinge. Clearly, a typical body may contain one or more hinge 
points, but a hinge may be associated with only two bodies. 

Hinge 1 is given special consideration. It is also a pair of 
points j but one of the pair is coincident with the reference 
point of body 1, and the other point of the pair is coincident 
with the inertial origin. Thus, motion "across the hinges" is 
used to represent system motion. The reference point of body 1 
is located with respect to the Inertial origin by an inertially 
referenced position vector. The attitude of the reference frame 
of body 1, with respect to the Inertial frame, is represented by 
three Euler angles. Thus, there are six position/attitude co- 
ordinates associated with hinge 1. 
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Each of the remaining hinges is considered in a manner somewhat 
similar to that of hinge 1. Referring to Figure 1-2, we note 
that there is an orthogonal reference frame attached to point p 
and another to point q. The triad of point p may have a ”natu- 
ral*’ (or undeformed) misalignment with respect to the triad of 
body point m plus additional misalignment due to ^lastic_^-_ 
formation. The same relationship is true concerning the points 
n and q . 

Now there are, associated with the hinge of points p and q, six 
relative position/attitude coordinates. Point q is located from 
point p with a p-frame referenced position vector. The attitude 
of the q-frame with respect to the p-frame is represented by 
three Euler rotations. Thus, if NH is the number of system hinges, 
then there are 6 x NH position coordinates to be used in conjunc- 
tion with modal displacement coordinates to define the system^ s 
position state. Let it be noted that only the time variable 
position coordinates of the 6 x NH set of candidates are consi- 
dered as state vector elements (the position coordinates whose 
rates are constrained to zero are not included; however, the po- 
sition coordinates themselves need not be zero) . 

The system of bodies generally has a number of so-called "sensor 
points." We define a sensor point to be a structural hard point, 
which has a right-handed orthogonal reference frame attached, 
that is used for a variety of purposes. A sensor point may be 
used to enable the program system to monitor the position, atti- 
tude, or the rates associated with a specific structural hard 
point. For example, a rate gyro, a star tracking device, or 
other motion/position sensing device is physically situated at a 
sensor point. Also, a sensor point is used as a point of appli- 
cation of a force or torque vector (see Figure 1-2). 

The system of bodies may contain built-in momentum wheels, some 
of which are constant speed wheels and others are variable speed. 
The variable speed momentum wheels are motor driven; the shaft 
torque results from a given control law. Each momentum wheel 
of the system must be associated with a sensor point because, 
for a general flexible body, the gyroscopic coupling is influ- 
enced by elastic motion. 

As is indicated in Figure I-l, the system may be in a non-topo- 
loglcal tree configuration. The methods employed in this develop- 
ment are such that c losed loop configurati ons Lgenerallv referr- 
red to as nontopological ) may be considered. If every body of 
the N-Body system is rigid, then of course there may be no closed 
loops, because such a system has an indeterminate "load path." 
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To accommodate closed loops, the system must contain sufficient 
structural flexibility (compliance), and therefore modal dis- 
placement coordinates, that the kinematic equations of intercon- 
nection constraint are algebraically consistent. 

The program development is such that none, several, or all bod- 
ies of the N-Body system may be flexible. The system may be 
"forced" by such environmental factors as gravity, gravity gra- 
dient, solar pressure, thermal gradient, and aerodynamic drag. 

The computer program system described herein falls into several 
categories of capability: (1) synthesis and time domain solution 

of the nonlinear differential equations of motion of the com- 
plete spacecraft system idealized as a collection of intercon- 
nected flexible (or rigid) bodies, (2) linearization of the 
governing equations by numerical means, (3) time domain solu- 
tions of the linearized equations that describe perturbation 
response of the complete spacecraft system about some predeter- 
mined (calculated or user -specif ied) nominal motion, and (4) gen- 
eral frequency domain stability analysis corresponding to the 
linearized spacecraft representation. 
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II. 


EQUATIONS OF STATE - TIME DOMAIN SIMULATION 


A. INTRODUCTORY DISCUSSION 


The state equations governing the dynamic response of a system 
of interconnected flexible bodies, that may be actively or pas- 
sively controlled and that may be "forced" by environmental fac- 
tors such as solar pressure, gravity gradient, aerodynamic drag, 
etc. are presented here in a concise summary form as; 

[II-l] {U}j = [m]'^ ({G}j + [b]J {A}j, 

IH-2] {C}j - {U}j , 

[II-3] {&} =S[B], {U} . 

[II-4] {h = f ({3}. { 3 }, (U, (b, {<5}) . 

subject to the constraint equations 

[II-5] S [b]. {U}. = {i}. 

J ^ ^ 

In Equations II-l through II-5 the index J ranges from 1 through 
the number of bodies of the system. Ecjuatlons II—l through 
II— 4 represent n first order, nonlinear, ordinary differential 
equations while Equation II-5 represents m additional conditions 
of kinematic constraint. Thus, the dimension of the state space 
for a given system of controlled bodies is (n-m) . That is, 
there are n— m state variables required to define the configura- 
tion at any instant of time t. 


State variables of the configuration space include ^solute ve- 
locities, {U}., modal displacements Position CQ^'dinaLag 

(both angula r^ and cartesian position ) IM, and additional vari- 
able^ {6} that we will subsequently refer to as control vari- 
ables; they are variables associated with the differential equa- 
tions that define a given control law. However, they may reflect 
any other auxiliary differential equations that are necessary 
to define the overall feedback system; for example, they may in- 
clude thermal equilibrium states or other state variables neces- 
sary to complete definition of a state dependent environment. 
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The right-hand sides of Equations II-l through II-4 are func- 
tionally dependent on all the state variables and time, although 
the relationships may be only termed implicit at this point* 

Let it suffice that, in a way of introduction, a description of 
the nature of the governing Equations II-l through II-5 be given 
here, and that more explicit development and discussion follow 
in subsequent chapters* 


The Equations of II-l represent the dynamic equilibrium equations 
for the typical jt/z body of the system* They are of the form 
shown whether the body is treated as rigid or flexible. They 

state, in effect, that a deformation dependent mass matrix [m],, 

c . J 

postmultiplied by a vector of relative accelerations {U}j, pro- 


duces a vector of inertial forces that is balanced by all other 
state and time dependent forces {G} . and interconnection con— 

T ^ 

straint forces, [b]. {X}* The vector {G}j includes inertial 


forces due to centrifugal and Coriolis acceleration, as well as 
elastic restoring forces, damping forces, control actuator for- 


ces, and so forth* 


The constraint forces 



{X} are necessary 


in order that the kinematic constraint equations (II-5) are sa- 
tisfied; elemet^„oi^the vector {X} are actually Lagrange multi- 
pliers, evaluated and used in the solution process. 


The Equations of II-2 simply represent a selection transforma- 

tion, because the vector of modal velocities ^ subvector 

of The Equations of II-3, used to develop {3>, represent 

a kinematical transformation, transforming nonholonomic veloc- 
ities to time derivatives of position coordinates. Finally,^ 
the Equations of II— 4 are auxiliary differential equations tnat 
^Tce. user defined and may be used to implement control dynamics 
and other feedback effects. 


The constraint Equations of II-5 are kinematic conditions of a 
form similar to those of Equation II-3. In either case, we have 
a velocity transformation. We might term Equation II-5 an ac- 
tive set of kinematic conditions and those of Equation II-3 a 
passive set. The active set is used to calculate m of the de- 
pendent elements of the {U}^ vectors in terms of the remaining 

independent elements and the prescribed velocities {a}, some of 

may be zero and some user— defined functions of time. Thus, 
the constraint equations are of a general form because nonholo- 
noraic, rheonomic conditions may be so represented* Given that 
the {U}. vectors satisfy the required conditions of Equation 

II-5, then the position rates, {0}, may be evaluated via the 
passive conditions of Equation II— 3, resulting in a kinematically 
consistent system* 
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[II-6] 


Note that there are m equations of constraint represented by 
II-5, There are also m Lagrange multipliers in the vector {A}. 
Most often, in studies of dynamic systems, the Lagrange multi- 
pliers and the dependent velocities and accelerations are en- 
tirely eliminated from the governing equations. Such is not 
the case in our development. We have chosen to involve Lagrange 
multipliers in our equations for two reasons; (1) we wish to 
monitor t he multipliers^ a s a function of system .motion, as they 
are interconpection_forces and tornuas. and (2) for purposes of 
numerical implementation it is convenient to calculate and use 
the {X} vector in Equation II-l. The Lagrange multipliers are 
calculated by differentiating Equation II-5 and combining that 
result with equation II-l giving 


-1 


=1 Z-r [b]j lm]~^ [b]j 


{a} - 

j 


{U}j + [b]j [m]"^ {G} 


Notice the following functional dependencies: 


[n-7] [b]j - f {Uj ) . 

[II-3] =■ f {Oj ) . 

thus 

[II-9] {p} = f |{p}, (c), {U}| , 

[II-IO] {bj - f |{u}j j , 

{p}, {p}, {c}, {Cl, {5}; t j , 
(U, {U}, {6}; t ), 

[11-13] [m]^ - f ( (Uj ), 

[11-14] [b]j - f ({p}, {p}, {O, {b ) , 


[II-ll] {6} - f 

[11-12] {G}j - f 


thus 


[11-15] {A} * f |{5}, {p}, {U}, (b, {0}, {«}; t 

[II-16] and {U} = f |{C), {p}, {U}, {|}, {^, {6 

where, in the above notation, we mean that the elements of the 
matrices/vectors on the left are functions of the elements of 
the vectors on the right. The chronology, of evaluations i ndicated 
is that which must be followed in the solution process. 
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The differential equations of motion for the system are therefore, 
of the general form: 


[11-17] y^ = f|yi, Y 2 . •• , t j, 

and the state vector and its time derivative arc arranged as fol- 
lows: 


{U>1 

{y} = 

{U>1 

{U>2 


(U}2 



■ 

• 


• 



^^>NB 

{Ui 


COi 

{Oz 


ibz 

m 


* 

* 


• 

• 




6i 


01 

• 

02 


02 

• 


• 

• 


• 

• 



^NB 




62 


62 
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• 

• 


• 

• 




. — M 
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with NB the total number of bodies of the system, N3 the total 
number of position coordinates necessary to orient the system 
and N6 the total number of auxiliary (control) differential 
equations required • 
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Now, given that the {y} vector is known (numerically) from pre- 
scribed initial conditions or from numerical integration of {y}, 
the primary task of the solution process is to numerically es- 
tablish the {y} vector. The {y} vector is numerically (step by 
step) integrated so as to produce an incremented {y} vector, thus 
a sequence of time point solutions. 

In way of summary, a narrative description of the steps (numeri- 
cal evaluations) necessary to produce {y} given {y}, follows. 


The matrices [B] , and ^.re kinematic coefficients that de- 

pend on position and modal displacement variables, and are eval- 
uated as the first step. 


Now, if available numerical techniques (also computer software 
and hardware) were absolutely accurate, we would be assured that 
the {U}, vectors, resulting from numerical integration of the 

. J ' 

{U}j vectors, would satisfy the constraint equation This 


is not the case, therefore the second step of the solution pro- 
cess is to calculate the dependent elements of the vectors 

by using Equation II-5. In fact, due to anticipating numerical 
inaccuracies, only the independent elements of the vectors 

are obtained by nximerical integration. There are only n-m **in- 
tegrators" involved in the solution process even though all of 


the elements of the vectors are numerically evaluated (by 

use of Equation II-l) ; we have good numerical resolution in the 


independent {U}, 
{X}. ^ 


elements due to using the Lagrange multipliers 


A kinematically consistent system results from satisfying Equa- 
tion II-5. The {U}. vectors may now be used with the selection 

and kinematic transformations as indicated by Equations II— 2 and 

II-3 to produce (numerically) all the modal velocities and 

position coordinate rates {3} completing the third step of the 
process. 


Sufficient calculation has been completed to this point to then 
evaluate^ the control variable rates as per Equation II-4, pro- 
ducing { 6 }. During the process of calculating the { 5 } vector, 
all of the required control actuator torques (or forces) are 
calculated, because sufficient numerical information is avail- 
able. All of the constituents of the torques/force vectors {G} 


are now available and therefore numeri- 

cally evaluated, (refer to the functional expressions of Equa- 
tions II-ll through 11-14), which completes the fourth step of 
the process. 


j’ 
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With reference to Equation II-6, we note that there is now suf- 
ficient ntuaerical information to evaluate {X}, which is then used 

in Equation II-l to calculate the Wy completing the fifth and 

final step of the process. 


It is noted in the above discussions that the solution process 
may be carried out through completion, providing the state vec- 
tor is numerically known. At any step of a simulation, the {y} 
vector is known, of course, as the result of numerical integra- 
tion. The initial state vector is another matter. It is diffi- 
cult, if not Impossible, for a user to prescribe vectors 

that are kinematically consistent with the conditions of Equa- 
tion II-5; also, the nonholonomlc velocities of when con- 

sidered as a complete set, are of a somewhat abstract nature. 

The user i.s in a much better posture to prescribe initial values 

of { 3 } and { 5 } (the initial velocities that are physically mean- 
ingful to him). Thus, to initiate the simulation (that is, to 
create an initial state vector from information the user is in 
a position to prescribe) some preliminary steps must be taken, 
as follows. 


The user must prescribe initial values of the { 5)^1 (3}» 

• J J 

{ 3 } and {6}^vectorsj also the^variable speed momentum wheel spin 
velocities 0. Now, in that {a} (the prescribed position rates), 
are explicitly dependent on time and are always available, the 
kinematic Equations II-3 and II-5 may be used together to estab- 
lish initial values of all {U}^. The question inevitably arises: 

are the number of equations represented by II-3 and II-5 suffi- 
cient to solve for the elements of the ? Let us consider 

the typical {U}, vector. We note that there are six reference 


frame velocities in each namely. 


X' y 2' 

There are also six relative velocities associated with each 
hinge. Now, if the system is a topological tree configuration, 
then the Equations of II-3 and II-5 comprise exactly the re- 
quired number of equations to establish the reference frame ve- 
locities; that is, there are as many hinge points as there are 
bodies and even if every body were rigid, the system would be 
determinate. In this case, the initial sets of six reference 
frame velocities are computed via Equations II-3 and II-5; the 


U) 


V, and w. 


prescribed initial { 5 } vectors and momentum wheel spin veloci- 
ties are simply placed in the appropriate 'CUlj vectors, and the 

initial state vector is thus defined. 
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In the event that the system is not a topological tree configu- 
ration, then there are more equations (II-3 and II-5) to be 
satisfied than there are reference frame velocities (or in other 
words, there are more hinges than bodies). In this case, ele- 
ments of the vectors must take on the responsibility of 

helping to satisfy the kinematic conditions. For each hinge 
in excess of the number of system bodies there must be at least 
six deformation modes, represented by ? coordinates, and they 
must be distributed throughout the system in such a way that 
the kinematic conditions of Equation II-5 are independent. 
Clearly then, w hen^ there are more hinges than bodies ( nontopo- 
logical tree) , one or more of the bodies must be flexibl e for 
the system to be determinate. Now, when the configuration is 
nontopological , the user will specify initial values for all of 
the 5, but he must acknowledge that they are not all independent 
and the dependent ones (automatically determined by the program) 
are calculated and replace the values that he has specified. 

From these considerations, we note that the initial state vector 
is created by the program from information that is user supplied 
and that is physically meaningful to him. His only concern, re- 
garding initial conditions, is: whether he has supplied an ade- 
quate description of system flexibility, in the event of a non- 
topological tree configuration, for the system's klnematical 
equations to be determinate. 


B. DERIVATION OF EQUATIONS OF DYNAMIC EQUILIBRIUM 


The differential equations of motion and auxiliary equations that 
characterize a physical system may take any one of several equi- 
valent forms. By equivalent form, we mean that the same physi- 
cal system can be characterized by more than one set of mathe- 
matical variables; in any case, the number of variables must be 
the same. For example, the motion equations for a rigid body 
might be derived by using Lagrange's equations (resulting in 
six second-order equations), or one might use the Newton-Euler 
equations where translational motion is represented by three 
second-order equations while rotational motion is represented 
by six first-order equations (three moment -moment urn equations 
and three attitude equations). In each case, there are 12 state 
variables • 
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There are a variety of alternative methods of analytical dynamics 
that one may select from to develop his final (programmable) 
equation format. A timely and valuable commentary accompanies 
the comprehensive comparative evaluation of these methods in a 
recent report by Likens’^. The basis for our development is ef- 
fectively included in his discussion. 


Our intent is not to highlight any particular method of analyt- 
ical dynamics as being superior to the others. Clearly, the 
methods are all equivalent providing that they are developed 
through completion without any compromising simplifications. 

The choice of method is made after considering the requirements 
associated with a particular problem or computer simulation pro- 
gram. Our development begins with a Lagrangian approach, then 
through algebraic manipulation we arrive at the format of Equa- 
tions II-l through II-5. 


[11-18] 


Lagrange equations for the general situation appear as 



9D 

■ 







X 


i 


n 


E 


‘ij 


q . + a . ^ * 0 

J It 


for n) 


for m) 


In these equations, T and V are system kinetic and potential en- 
ergies, respectively, and D is the Rayleigh dissipation function 
(accounting for internal damping) • The generalized constraint 

forces^ ^ j a^^ augment the generalized forces Qj (that arise 

due to the action of external factors) and are necessary in order 
that the additional conditions of constraint (the second set of 
Equation 11-13 be satisfied. The form of the Equations 11-18 
is complete and general, in that they include unconservative 
forces (explicitly time dependent Qj and disslpitive forces 

3d/ 3qj and the auxiliary constraint equations (the second set 

of Equation 11-13) are in an all encompassing form, because 


*Llkens, P, W. , "Analytical Dynamics and Nonrlgid Spacecraft Sim- 
ulation," Technical Report 32-1593, Jet Propulsion Laboratory, 
Pasadena, California, July 15, 1974. 
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holonomic conditions may be so represented. The coefficients 
(a.., 2, n; t) may depend explicitly on the time (t), 

thus the constraint conditions as shown account for both rheo- 
nomic and scleronomic situations. 

In the equations, n is the number of generalized coordinates in- 
volved in the representation and m is the number of auxiliary 
conditions of constraint. Note that, although the are gener- 
alized coordinates (as they must be for the Lagrangian formula- 
tion) they are independent only in the isolated case when m=«0, 
or when there are no aiaxiliary constraint conditions. The wri- 
ter has observed that some engineers share a misconception on 
this point, thinking that if the variables q^ are not indepen- 
dent then they are not generalized coordinates. In view of the 
m constraint equations, we simply have a set of generalized co- 
ordinates that are not independent. 

In cases where all of the constraint equations are holonomic, 
it is theoretically possible to eliminate m of the q^ in terms 

of the remaining n-ra. However, if any of the constraint condi- 
tions are nonholonomic , a Lagrange multiplier must be used 

in conjunction with that equation. Lagrange multipliers may, 
of course, be used for either holonomic or nonholonomic con- 
straints • 

In that the simulation program includes mathematical representa- 
tion of active or passive control for elements of the spacecraft 
system, there are state equations involving control variables 
that are additional to 11-13. The manner in which the additional 
control equations enter into the composite system state equations 
is the same whether we are talking about the form given by Equa- 
tion II-l or that of Equation 11-13. The control system state 
variables retain their identity in either case although the con- 
trol forces/ torques necessary to **close the loop” are trans- 
formed differently. In the case of Lagrange's equations, the 
control torques contribute to the generalized forces whereas 

in the case of the summary Equations II-l, they contribute to 
elements of {G} and may be interpreted to be ordinary forces or 
torques, acting at a structural hard point (or at a sensor point). 
Thus we will postpone further discussion of the control system 
until later, concentrating on the "mainline” motion equations 
until such a point when we can clearly indicate control system 
coupling. 
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In order to "solve" Lagrange’s equations of motion, one must first 
define the explicit form of the kinetic and potential energy func- 
tions, the dissipation function D, and he must also define the 
form of the transformation relating ordinary cartesian position 
coordinates (positioning the typical system particle or element) 
to the generalized coordinates the form of the transformation 

is necessary to be able to express generalized forces Qj in terms 

of external ordinary forces. Having defined the form of the en- 
ergy functions and coordinate transformation, one merely per- 
forms the indicated differentiations (11—13). He has not yet 
solved the motion equations but has only explicitly defined a 
system of ordinary second— order differential equations, which 
in many cases are nonlinear, and which require solution using 
numerical integration techniques. 

With numerical implementation and digital programming in mind, 
we wish to recast the form of the ordinary differential equa- 
tions. First of all, we would like for them to result in canon- 
ical first order form (the highest time derivatives appear un- 
coupled on the left hand side). Also, we would like to group 
complicated combinations of generalized velocities and displace- 
ments so that we may replace such groups with new variable names. 
The new variables we refer to have been called "quasi-coordi- 
nates" in the literature. This will simplify the required com- 
puter programming and minimize arithematic computation. Also, 
it helps considerably in organizing the numerical algorithms 
necessary to evaluate the left hand side of the state equations. 
Thus, recasting the form of the governing equations is suffi- 
ciently justified. 

We begin the recasting process by defining the forms of kinetic 
and potential energy, and the required transformation. First 
let us note that bodies of the system of flexible bodies are 
tentatively treated as though they were completely independent, 
one of the other. The influence of any one body on another is 
accounted for through the additional constraint conditions and 
the Lagrange multipliers. Thus, if we express kinetic and po- 
tential energies for the typical body and apply Lagrange’s equa- 
tions to it, the ordinary differential equations pertaining to 
it are simply a subset of Equation II— 13j and we will have ac- 
counted for the total system through the representative form of 
the typical body. 
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The generalized coordinates chosen to represent the configuration 
of the typical body Include three Euler angles to Indicate atti- 
tude of the body fixed axis system relative to an inertial frame, 
three projections (components) of the position vector from the 
origin of the inertial frame to the origin of the body fixed ref- 
erence system, onto the inertial axes, and N elastic displace- 
ment coordinates. We note that the origin of the body fixed axis 
system needn't necessarily coincide with the body's mass center. 
Also, the elastic displacement coordinates may be measurements 
of displacement at a discrete set of points on the body or they 
may be coordinates associated with normal vibration modes. In 
either case, they represent displacements measured in the body 

th 

axis system. For the r flexible body, we tabulate its gener- 
alized coordinates as; 


(q,} 



Attitude 
Euler Angles 


Body's Reference 
Point Position 
Coordinates 


Elastic Displace- 
ment Coordinates 


Now, there exists a transformation that relates a set of nonhol- 
onomic velocities to the generalized velocities that is exten- 
sively used in recasting the equations. The transformation 
appears as follows: 
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where in Equation II-19 the vector of nonholonomic velocities 
{U} contains the three projections (u>^, u)^) of the angular 

velocity vector U onto the body fixed axes (oT is the angular ve- 
locity of the body reference frame), the three projections of 
the reference point translational velocity (u, v, w) onto the 

body fixed axes and the displacement rates The elements 

of the transformation (i, j®l, 2, 3) are direction cosines; 

the subtnatrix [y] is an orthonormal rotation transformation re- 
lating the attitude of the body fixed axis system to the iner- 
tial frame. The submatrix [tt] is also a rotation transforma- 
tion; however, it is not orthonormal because it relates vector 
components based on an orthogonal basis to those of a skew (non— 
orthogonal) basis; namely the axes about which Euler rotations 
are measured. 

In short, we write 

[11-20] {U} = [0] {q}. 

Clearly the elements of [0] are functions of the three Euler 
angles. There are 12 possible sets of Euler angles. Any one 
set is valid for use in subsequent development; the resulting 
equation form is independent of selection from the 12 sets of 
angles. 
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Elements of the transformation [3] 'stsiy be explicitly defined in 
terms of three of the generalized coordinates (the Euler angles). 

The kinetic energy expression for the rth body is most easily 
expressed (initially) in terms of the nonholonomic velocities^ 
{U}. Having done this, [3] is used to replace {Ul with [3] (q). 
The kinetic energy is then expressed completely in terms of gen- 
eralized displacements and velocities (the form necessary for 
applying Equation 11-18) . 

Kinetic energy for the typical body is 

[11-21] T ^ h f v*vadV 

•'v 


where v is the velocity field, a is mass density, and where in- 
tegration is carried out over the volume V of the body. 

Tae inertial position of any point p of the body is (See Figure 
II-l.) 


[11-22] r * + -po + n 



(R, the origin of the body axis system), positions the point 
p' (which coincides with p in the undeformed configuration) from 

point R, and where IT (x, y, z, t) is a measure of elastic dis- 
placement. 


The vectors pq and n are referenced to the body axis system, 
thus 


[11-23] ?o 


1^1 J k 


X 

y 


and 

[11-24] "n (x, y, z, 




*yk 



the elastic displacement n is represented as the superposition 
of a finite number of single valued space functions cj>^. 
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The velocity field v is obtained as 
[II-251 V - ^ + ¥ X ( po + n) + 

k=l 


with 


V = ‘'"r. 


dt 


The velocity of the reference point R may be expressed in terms 
of components referenced to either the inertial frame or the 
body frame, that is 


[11-26] "^R - J K 


also 
V 


K -II T il 


w 


The unit vectors {i, j, k}, {I, J» K} are related through the 
rotation transformation ly] follows that 


[11-27] 


u 


w 


= [y] 


At this point, let us introduce the repeated index summation con- 
vention to be concise. With this convention, when any two fac- 
tors of a term have the same index, summation over the range of 

that index is implied and thej^sign is deleted. For example, 
the third terra on the right of Equation 11-25 is 

\ K 

and represents 


**k ^k • 


Now, if we substitute H-25 into H-21, the kinetic energy is 
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[11-23] T 




+ n? X (po + n)] • [w 


^ \ 

+ 2 V • [u" X Cpo + n”)] + 2 V • ■?■ i 


k ^k 


+ 2 [uJ X Cpo + "n)j * I 

or, integrating term by term over V, 
[11-29] T * ^ m j^u V wj {u v w} 


+ ^ ju) 0) CJ I 

[ X y z J 


®jk ^k 


J -J -J 

XX xy . xz 

-J J -J 

yx yy yz 

-J -J J 

zx zy zz 


j^u V wj 


j^u V 


0 s -s 

z y 

-SOS 
z x 

S -S 0 

y X 


xk 


X 


yk 


zk 


+ |w W 0) 1 
L X y zj 


xk 


yic 


zk 


where we have used 


[11-30] 


"■/ 

•'tT 


QdV 


XX 






[11-31] 


J + 2(b , + b .) ^. + (c . k 

xxo yyj zzj J yjy*^ 


X (P'0 + “iD] 


odV 

^zjzk^ ^k 
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with 


[11-32] b 


yyj 


zzj 


and 


/ 

I 

'I: 


y 4>yj odv. 


z , adV, 
Zj 


c . , * / iji , , adV 

yjzk / ’^yj ’^zk 


Also, we have used 


[11-33] a 


xk 


tII-34) 


L 

li 


‘*>xk 


^xj “^xk '^yj “^yk 


^ ^k) 


adV 


and 

[11-35] S 




(" *yj '3 

^yzk * '’zyk ^ ^^^yjak *^zjyk j ^j’ 
and also, 

■^xyo ^|*’xyj ^ Nxj) ^3 ^ ^x3yk ^3 ^k‘ 

All other quantities involved in Equation 11-29 are obtained by 
cyclic permutation of the indexes x, y, and z. Finally, as the 
kinetic energy is of quadratic fom in the elements of {U}, we 
may express it as a triple matrix product 


)*xk ^3 S ) *yk 


[11-36] d 


xk 


/ 
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[11-38] T = % w [m] {U} 
with 

[11-39] [m] = 


J -J -J 
XX xy xz 


-S s ' d , 
z y 1 xi 

"x 2 

• • 

^xN 

J -J 

yy yz 

S 

z 

-Sid, 
X 1 yl 

y 2 

• • 

**yN 

J 

zz 

-s 

y 

^x 1 

‘^za 

• • 

^zN 


ni 

1 ®xl 


• • 

®xN 



m 1 ^yi 

a 15 

• • 

®yN 



* a , 

a 

• • 

a „ 



m 1 

z 2 


zN 


— - 

1 

— — 

. > 


(Synunetric) 


|eii 

1 

^12 

• • 

® 1 N 



1 

1 

«22 

• • 

® 2 N 



1 

1 


• • 

®NN 






—* 


or in short, 
[11-40] [m] 


i W 1 

S a 

T • r 

d I a 1 a 


Using Equations 11-40, 11-19, and 11-38 gives 
T 


[11-41] T * ^ 


Clearly, the elements of [m] depend on only the the elements 

of [ 8 ] depend on the Euler angles and therefore kinetic energy is 
a function of generalized velocities and the generalized coordi- 
nates themselves j thus, the functional notation 


X - T ^qi, q2» 

is applicable; 
important role 


•• q^; qi. 92, •• <ln) 

terms such as ^T^'Sqj will come about and play an 
in the simulation. 
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To continue it is necessary to express the potential energy V 
and dissipation function D. Let us assume that the elastic 
strain energy can be written as a positive-definite quadratic 
form in the elastic displacement coordinates, or 

[11-42] V - Js [k] U); 

the syncnetric matrix [k] is developed by standard finite ele- 
ment techniques such as those embodied in NASTRAN, In the event 
is a set of normal modal coordinates, then [k] is diagonal 

th 

with the j diagonal element appearing as 

rii-43] 

with w being the natural frequency. Of course, normaliza- 
tion of the eigenvectors (mode shapes) is assxnued such that the 
generalized mass for the j vibration mode is unity. 

Now, since 

[11-44] ^ [0|0!I^] {q} 

= [S^l {q} 
it follows that 

[11-45] V = *5 jjlJ [S^]^ [k] [S^] {q}. 

Similarly, D is written as 

[11-46] D = is [qj [S^]^ [C] [S^]{q}, 

the matrix [C] being equivalent viscous damping for the struc- 
ture; it is also developed using standard finite element tech- 
niques. 

Let us now refer back to Lagrange^s Equations (11-13), and re- 
express them in matrix format 

[0 .]^ [m] [0] U)! 

. ) * 

+ ij I |^<lj [^] [®] j^qj [0] [<l}|+[a]'^{l} 


+ {Q} + is 
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and 


[11-48] [a]{q} - • 

What is meant by [B J and [m J is the partial derivative of 

> J *3 

every element of [B] and [m] with respect to the j generalized 
coordinate. 

Let us now define the ordinary momenta 
[11-49] {p> - [m] [8] {q} 

- [m] {U} . 

Also, since {U} = [8] {q} 

[11-50] it follows that {q} = [8] ^ {U}. 

Using Equations 11-49, 11-50, 11-47, and 11-48, we may write 


[11-51] 


(p) - [s^] {q) + [c] [Sj] (4) 


4 .) 


and 


[6]"^^ {Q} + [3r^^^j[qJ {p}j - [B]"^ {p)j 

|LuJ [m^j] {U}| + [8]“^’^ [a]’^ {X}, 


+ ^ [ 8 ] 


[11-52] [a] [8]"^ {U} - {-a^}. 

Several observations can be made on studying Equations 11—51 and 
11-52: 

First of all, recall the form of [8] and [S^] (Equations 11-19 

and 11-44). It is clear from these forms that 

— IT T T 

[11-53] [8] [S^]" = [S^] 

[11-54] and that [S^j {q> - 

[11-55] and [S^] {q> -{b. 

Also, since the elements of [m] depend only on the first six 

elements of ) LuJ [m j] {U}(are null, thus 

[11-56] [8]"^'^) LUJ [n>^j] {U}| = |luJ Im^.] {U} j . 
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-iT 

Further, we note that the matrix [3] transforms the general- 
ized forces {q} to forces ^‘acting in the quasi-coordinates,” or 
let us call 

[11-57] {G^^} = [0]“^^ {Q}. 

ex 

thus contains ordinary forces and moments due to external 

sources and corresponds to time derivatives of the ordinary mo- 
menta. 


Because the transformation [&] depends only on the Euler angles, 
it follows that only the first six elements of the column 

[0]‘^'^(| [qj {p}j - tpij 

are non-zero, and one finds after considerable algebraic manip- 
ulation that this column may be reexpressed as 

in] {p} 


or 

[11-58] W {p} 


0 u) -uj 1 0 w -V 1 


p(w ) 

z y 1 I 



-w 0 0) ! -w 0 u 

Z X 1 1 


P(“y) 

oj -0) 0 , V -u 0 


P(w^) 

y X 1 .1 



1 0 m -oj 1 


p(u) 

1 z y ' 



1 ) 
, -0) 0 i 


p(v) 

\ Z X 



o 

3 

1 

3 


p(w) 

» y X 



- -i “-1- 

I 1 

1 1 


■p^l” 

• 

1 1 
1 1 


« 

.P(Sy)_ 


With these observations and definitions, the Equations 11-51 and 
11-52 may be reexpressed as 


[H-59] (p) - 

+ ^ + Ib]’^ {A}, 

[11-60] and [b]{U> - {a} 


(U + [jj]{p} 
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where we have used 


[H-61] [b] = [a] [3] ^ 

and 

[11-62] 

Notice that the constraint equations (11-60) are now expressed In 
terms of the nonholonomic velocities {u}; the coefficients [b] 
are obtained directly from relatively simple, vectorial expres- 
sions of kinematic constraint. The same [b] coefficients are 
transposed and used to multiply {X}, producing constraint forces/ 
torques corresponding to the ordinary momenta. 


[11-63] 


[11-64] 

to be used in conjunction with system kinematic constraint equa- 
tions 

[II-65]2[b]j. {U}^ = M 
r 

which is the same form as that given by Equations II-l and II-5. 

The last three terms of {G} given in Equation 11-63 are inertial 
forces that involve velocities and displacements of the body. 

The matrix [m] is an instantaneous inertia matrix, depending on 
instantaneous values of the deformation coordinates The 

centrifugal and Coriolis effects are completely accounted for 
within the framework of the assumed velocity field (given by 
Equation 11-25) . These effects would not be accounted for if 
we neglected "tangential" velocity due to elastic displacement; 

that is, if we assumed that ju>xn|<<lwxpol- In this case, the 
inertia would be constant, independent of * 


If we now define the {G} vector to be 


{G} 


{G } - 
ex 


{ 5 } - 


(b + [fi] [">] tu> 


+ IjJ) [m j]{U} j- [m]{U} 

it follows that we may write dynamic equilibrium equations for 

th 


the typical r body as 

{U}^ = + [b]^ {A}j 
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An accurate definition of the dynamic equilibrium equations 
clearly hinges on a complete and accurate definition of the con- 
stituents of the vector, which includes the inertia matrix 

[m]^. Also, the kinematic coefficients [b]^ must be developed 

in an exact fashion. Kinematics and a more explicit development 
of {G} are given in subsequent sections. 
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KINEI-IATICS AND SYSTEl'I TOPOLOGY 


From a Lagranglan formulation all of the generaln^zed forces, not 
derivable from a potential function, ordinarily appear as {Q} on 
the right side of Lagrange equations of motion. We have accounted 
for internal damping forces with the use of Rayleigh’s dissipation 
function D and for generalized constraint forces through use of 
Lagrange’s multipliers. 

Thus, the generalized forces that remain to deal with include those 
due to external factors such as aerodynamic drag, solar pressure, 
and other commonly encountered environmental loadings. 

We also intend to treat control forces (servodrive torques, reac- 
tion jets, etc.) as though they were external. They are not ex- 
plicitly external, of course, because they depend on time through 
position and rate errors that are functions of elements of the 
state vector and on control system state variables that arise from a 
given control law. 

Let us assume that there is a finite number of points on the typ- 
ical body where a force vector (or torque) is known to act. Each 
of these force/ torque vectors contributes to the generalized 
forces {Q}. The generalized forces are calculated by expressing 
the virtual work of the external ordinary forces in terms of vir- 
tual displacements of the points of force application. The trans- 
formation relating ordinary coordinates to generalized coordinates 
is then used to define the explicit form of the generalized forces. 

For example, suppose that a force f^ and torque T^ act at point p 

of the typical body. Their virtual work is 

661 6W - ? • 67 + T • 5? . 

P p P P 


Notice that we treated the virtual rotation 60^ as a vector quantity. 

This is valid, even though a general rotation is not a vector quan- 
tity, for the virtual rotation is infinitesimal and therefore is 
a vector. Further, because virtual displacements are Infinitesimal, 


we may express 


6r and 60 in terms of virtual displacements of 
P P 

the quasi-coordinates; that is 



[11-67] 6r 


[1 j ij ( 


and 


[11-68] 


[i j kj( 


6ri 

r 

6r 2 
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6x2 


*^xj 

(Xp, y 

*yj 

(Xp. y 
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CD 
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o . (x , y , 
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^ A (x , y , 
zj p’ ■'p’ 
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Z ) 
p 






where (6ri, &V 2 » <Sr 3 ) are components of virtual displacement of 
the body's reference point R» components of 

virtual rotation of the body axis system, and (o , _ a ,, a 

th L yy zy 

are components of the j space function representing elastic 
rotation at point p (modal slopes, for example). 


Now, let us assume that the force and torque vectors (T and T ) 

P P 

are referenced to the body axis system, thus they may be written 
as 


[11-69] f 


J "J 


xp 


and 


[11-70] Tp = [i J 


yp 


zp 


xp 


yp 


zp 
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We note that virtual displacements of the quasi-coordinates are 
related to virtual generalized displacements by the same trans- 
formation that relates nonholonomic velocities to generalized ve- 
locities (See 11-19). It follows that the virtual work due to 


f and T may be written as 
P P 



The virtual work is also expressed 
6W - [6qJ {Q}, 


and because 6q^ is arbitrary and independent (it is treated as 

though independent in the face of Lagrange multipliers and con- 
straint equations) it follows that 


[11-72] {Q} 


Ibp]^ 


fti- 


The Equations 11-71 or 11-72 have a noteworthy geometrical in- 
terpretation. Notice that the first three lines of [b^] | 

are components of the resultant torque vector T^+(po+n) x f^, 
acting at the body’s reference point R, The second three lines are 
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components of the resultant force vector f^, while the j 


th 


line 


(j>6) corresponds to the standard procedure (of structural dynalI^- 
icists) to calculate it is usually expressed, general- 
ized forces acting in deformation modes are 


{Q} - [^] {f}. 


Also, recalling the form of [3], (Equation 11-19), we note that 
T 

[tt] resolves the resultant torque vector (about orthogonal body 
axes) to components about skew axes about which Euler rotations 

T 

are measured while [yJ resolves the resultant force vector (about 
orthogonal body axes) to components along the inertial axes. Further, 
we notice that [b^] is a matrix of coefficients that relates the 

velocity of any point p to the vector {U}. This gives us some ad- 
ditional insight as to why the same coefficients that are used in 
the kinematic constraint equations (11-60) are used (in transposed 
form) to multiply {X} producing resultant constraint forces. 


Thus, we have pointed out the remarkable duality of purpose asso- 
ciated with [b] type coefficients. They are initially expressed 
by writing simple kinematic velocity relationships. The coeffi- 
T 

cients [b] are then used to transform discrete ordinary forces 
and torques to equivalent forces and torques acting through the 
body’s reference point R. The matrix [3], which is also a velocity 
transformation, is transposed to produce the transformation to 
generalized forces (should they be desired). 


For our ordinary momenta equations we simply wish to express {G } 

ex 

which (following Equation 11-57) is given by 


[11-73] {G } 
ex p 


{Q} 

l&f [bp]^ 


Ct 




This {G } given by 11-73 reflects only the contribution of the 

force/torque acting at a single point p. The total must be 

obtained by summing over all the points of the body where forces 
and torques act, or 
NP 


[H-74] {G^^) 


E 

1=1 


Ib„ ]■ 
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Kinematic coefficients [b ] such as those of the previous example, 

P 

will be required throughout in our formulation of the state equa- 
tions, They are used to synthesize the constraint equations, to 
produce {G}, and they are even involved in the velocity transfor- 
mation of I 1-3, It is therefore advantageous for us to think of 
a ”bank” or collection of all the required kinematic coefficients 
to be put together in a semiautomatic fashion by using input 
specifications to the digital program. 

1. Sensor Point Kinematics — Force/Torque Transformations 


Consider the typical structural hard point s (See Figure II-2) • 
Let us assume a right-handed triad is fixed to point s and that 

the elements of the triad are unit vectors labeled T, m, and n. 
Now body n (which has point s on it) also has a right-handed 
triad fixed to point n. Suppose that, even when body n is in an 
undeformed state, the s-triad is misaligned with respect to the 
n-triad. When the body deforms there may be further angular mis- 
alignment between the two triads. Thus, the relationship linking 
the two sets of unit vectors is 


[ri-75] 

■“ " 



I 


T 


’m 


T 


n 


k 


with [ R and [ >iR ] being orthonormal rotation transformations, 
s s s n 

the first relating the "naturally" misaligned triads via constant 
Euler rotations and the second accounting for additional rotation 
due to the body’s deformation at point s. 


The structural deformation at point s is assumed to be sufficiently 

small that the Euler rotations associated with [ ] may be eval- 

s n 

uated through use of 


[11-76] 


9l 

02 

03 


[O^] {C>, 


where [o ] is a (3xN) matrix of modal rotation amplitudes (each 

. .s ’ 

of the N columns corresponds fo a deformation mode) at point s. 
Let us consisely denote the triads associated with points n and s 


by {e } and {e } respectively. Then we may express the relation- 
n a 

ship linking the two sets of unit vectors as 


[11-77] 


{¥ } = 


R 

s n 


{e }. 
n 
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There is a requirement for expressing the absolute velocity of a 
typical s-point and the angular velocity of the typical s-triad, 
in subsequent kinematic development, in terms of velocity states 
of a given body. Let us think of a six long vector (column) of 
velocity components (three rotational and three translational) 

that are projections of m and v onto the s- triad axes* It is 

s s 

related to the vector for the body by the transformation 



with [h ] and [o ] representing matrices of displacement and ro- 
tation amplitudes, respectively, and with being an anti- 

symmetric matrix accounting for a vector cross product, or 


;(n) 

ns 


-(2 +T1 ) 

s zs 

y +n 

ys 


2 +n 

s zs 




X +n 

s xs 


-(x +n ) I 
S xs I 


The superscripts used in Equations 11-73 and 11-79 are used to 
indicate the frame to which the velocity components are referenced. 
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Kinematic coefficients such as those of Equation 11-78 are gen- 
erated for each so-called sensor point of the system of bodies. 
They are used by the simulation program to produce contributions 
to from given force/ torque components in the manner indi- 

cated by Equation 11-74. 


2. Hinge Point Kinematics 


Kinematics associated with hinges follows a line of development 
somewhat similar to that of sensor points. Consider the points 
p and q (refer to Figure II-2) to be two structural hard points 
associated with a given hinge. All necessary kinematics infor- 
mation pertinent to the hinge is obtained through expressing the 
velocity of point q relative to point p and in expressing the 
relative angular velocity between the q and p frames. It is conven- 
ient that the angular velocity components are projections onto 
skew axes (Euler angle rates) and that translational velocity com- 
ponents are projections onto the axes of the p triad. Let us as- 
semble the six relative velocity components into a column matrix 
as 


[11-30] {3), 


{9} 

{A} 


with being the three relative Euler angle rates and being 


the three relative translational velocity components all pertain- 
ing to the k hinge. Now tne column of relative velocities may 
be expressed as 


[11-31] 

with 


^ «'n 


[11-82] [bp] - 

[q“p] [p". 

11 i - [A] M pp] 1 


[p\] 

r 1 ' M ["■’] 

J 

[11-33] [bq] - 



• 


[pN] 

jP'',] [q®n] 1 [p%] [q®n] [\] J 
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In Equations 11-82 and 11-83 the rotation transformations [^R ] 

^ pm 

and developed to include the effects of structural de- 

formation in the sense indicated in Equation 11-75; the rotation 
transformations [r] ^ and developed in standard fashion 

using the three Euler rotations 


NOTE: 

Hinge labels are circled; 
body labels are not circled. 



J 


Figure II- 3 Topology of a Typical System 


For purposes of further discussion, consider the system of bodies 
of Figure II-3. Topology of the system is simply indicated by an 
integer array we call ITOPOL, which is as follows: 

12345678 ^ Hinge number 

■Body (n) relative to 

•Body (m) 


The [ITOPOL] array, which is actual input to the simulation pro- 
gram, is used to define system topology as indicated* Now, with 
reference to the example shown in Figure XX— 3 and the correspond- 
ing (ITOPOL) array, let us indicate the form of the velocity 
transformation. We may write 


[ITOPOL] - 

12435677 


03263152 
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where ] and [b^ ] are matrices as defined in Equations 

11-81 and 11-83 (with i=Hinge number and j=Body number). The ve- 
locity transformations of Equation 11-84 represent the "bank" of 
all hinge kinematics coefficients previously mentioned, and pro- 
duces every possible velocity component pertinent to hinges. Re- 
ferring to the basic system equations II-3 and II-5, we note that 
selected lines, or equations, from the bank (11-84) are taken to 
represent constraint equations or position coordinate rate equa- 
tions. The [B]j and [bj^ coefficients of Equations II-3 and II-5 

are simply subpartitions extracted from Equation 11-84. 
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To implement calculation of Lagrange multipliers (refer to 
Equation II-6) it is necessary to develop time derivatives of [b]j 

coefficients. In a manner similar to above, where all 

efficients are extracted from the complete collections, the [b]j 
matrices come from a collection of matrices whose members are 


[b ] and [b 

‘ii.j 


] 

i.j 


which are developed In Appendix C. 
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D. DEVELOPMENT OF THE {G}^ FORCE VECTOR 

The equations of dynamic equilibrium for the j body of the sys- 
tems are given in an earlier section as Equations II-l. As was 
noted there, the right-hand side includes a so-called {G}j vector, 

which accounts for all state dependent forces except for those 
of interconnection constraint. Earlier in Chapter II (Equation 
11-63), the {G} vector is presented in a somewhat more developed 
form. ^ 


The purpose of this section is to provide more explicit develop- 
ment of the elements contributing to {G}^ , Let us account for 

all contributions in the following expression (we omit the J sub- 

script, understanding that we are dealing with the typical, or j 
body) ; 



"o" 


"b’ 

III-85] {G} - {G^^} - 
/ 

k 

{U - 

c 


+ 


I L?J - ‘"1 * '“gg’ 


The first term {G } has already been discussed in the previous 
ex 

section (See Equation 11-74), but we note here that the ordinary 
force/torque components that produce thought of as a 

miscellaneous force vector. Its presence provides the program 
user latitude to include a variety of additional effects. Clearly, 
it is the Implement through which control forces/ torques are 
”fed back” to the dynamic system. 


The second and third terms of Equation 11-85 have been previously 
introduced. There is no implicit restriction on the stiffness and 
damping matrices [k] and [C], nor is there a restriction on defi- 
nition of the {?} coordinates; they will likely be coordinates 
associated with orthonormal vibration modes in the majority of 
cases. However, they may be physical (ordinary-discrete) displace- 
ment coordinates as well. In the latter case, the [k] and [C] 
matrices are generally coupled. 
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The last two terms of Equation 11-85 are included to account for 
momentum wheel coupling and gravity effects respectively. The 
treatment given to built-in momentum wheels is such that, in ad- 
dition to producing a contribution to {G}, there is also a re- 
quired extension to the form of the [m]j matrices. This is be- 
cause momentum wheels are ■vn&Tti-a.Vly coupled. Thus, there is 
sufficient requirement for a dedicated development concerning mo- 
mentum wheels. The following two sections deal exclusively with 
momentum wheel and gravity effects, respectively. 

The remaining terms contributing to {G} are basic inertial effects 
and involve the matrices [mj , I™ » and Im] . With reference to 

Equation 11-39, the form of [m] is given corresponding to the 
case where one has single valued space functions cjij^ available to 

him. Ordinarily, one does not have access to such a description 
of the structure's deformation modes, due to the structural com- 
plexity of typical spacecraft. The analyst should always be able 
to obtain, .as data, matrices of modal amplitude ratios ("mode 
shapes") and the corresponding structural mass matrix (generated 
by use of finite element techniques) . To accommodate data based 
on the more practical definition of structural characteristics, 
it is necessary to recast the inertia matrices [mJ in a similar 
but more general format. The generality of the development of 
Section II. B is not compromised by extending the form of the in- 
ertia matrix. The extended, or more general, inertia matrix is 
developed in Appendix A, but here, for purposes of developing 
inertial contributions to the {G} vector, let us accept the re- 
sulting form; and present the kinetic energy expression as 

[11-86] T - *5 j_uj ^[m^] + [mjj j 

with the repeated index summation convention implied, and with 
[m ] of the form 

[11-87] [m^] » 


that is it. is just like the [m] given by Equation 11-39 except 
it is constant, independent of deformation. The constant iner- 
tia matrix [nip]» as given by Equation 11-87, Is always of the 

form shown regardless of the choice of "modal" coltnnns. The form 
of the matrices [m;] and [m 2 ] is Puch as to accommodate the gen- 
eral situation; that is, their definition includes Inertial inte- 
grals as defined for a continuous system, (Equations 11-30 through 
Xl-37), or as defined by structural mass matrices that are called 
"lumped" or "consistent." 


-S 


S I m 
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The inertia tnatrix associated with Cj is 



and the one associated with is 

[11-89] “ 


jk 

Now, for N deformation modes associated with a given body, it is 
understood that the range of the indices j and k is N, thus the 
coefficients (Cn),, , (Ci 2 ) 4 i,i * ‘ (C ),, are stored as 9 (NxN) 

arrays of inertial integrals while (bi)^, (b 2 )j» •• 

(ai)j, ( 0 ^ 2 ) j» stored as a (6xN) array and a (9x^0 

array respectively. Thus, from a programming standpoint, we 
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note that there are 9N^ + 15N storage locations required to ac- 
commodate the inertial integrals necessary to account for the 
deformation dependent mass matrix. Of course, if a particular 
body is rigid (N=0) then only the first ( 6 x 6 ) diagonal partition 

of [m 1 is used, 
o 


l^en the body is flexible (N>0) then the inertia matrix is cal- 
culated from deformation states (Cj)and inertia Integrals in the 

manner indicated by Equation 11-86; the redundant operations due 
to symmetry and null operations are avoided in the digital code. 


Having an instantaneous numerical evaluation of the Inertia ma- 
trix, the term [Q] [ra] {U} is calculated and added to (G), con- 
sistent with the expression of Equation 11-58. 

It is now possible to express explicitly, the combination of the 
remaining two inertial force vectors in terms of the inertial 
integrals given in Equations 11-88 and 11-89. For purposes of 
further development, let us define the combination as 


[11-90] {G^} - |[UJ [m {U} j - [m] {U>. 


Thus, the first element of corresponding to is 


[11-91] (G^)i 


-2<o^ (bl)j + “z 


-u (ai)j - V (a 2 )j - w ( 03 )^ 




the second element, corresponding to is 
[11-92] (G^>2 = I (bi+)j -2Wy C^2)j + “2 

-u (Q4 )j - "V (“5)j - W («6)j 

-2Wy (C22)jij + “2 + <C23)jji] 

the third element, corresponding to w is 

z 
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[11-93] (G ^)3 - jo)^ (b 5 )j + (bg)j - 2w^(b3)j, 

-u (“?)j -V (a8)j “W (oi9)j 

■^Vjk ^ “x h 

+ “y + (C23)j,i] -2u^ <C33)^j j S 

the fourth element, corresponding to u is 

[11-94] (G >4 = -L (aj) + u ( 04 ) + to (ay) 4 ! 4. , 

^ jx j y j z jij 

the fifth element, corresponding to v is 
[11-95] (G^)s = ~|“x ^“ 2 )j + t^y ^“s)j + “2 

and the sixth element, corresponding to w is 


[11-96] (G )g = - { (u ( 013 ), + u (ag). + “ («9) 


j 


j 


'j 


1 1 . . 
( ^ 


Finally, for the element k-f6, corresponding to an inertial force 
acting in the C^coordinate we have 

[11-97] ^ ^2 [(Cii)^j Cj + (bi)j^] 

+ + <b2>j^] 

+ co2 [(C33)j,j Cj + (ba)^] 


- { [(Ci2)^j + (Ci2)^^] C, + (b4> 


X y 


-WO) 
X z 


0) OJ 

y z 


"j 


[(Ci3)kj ^ 


[(C23>kj (C23)jk^ 

\ I \ — \ 1 f 


+ 0 )^ ^ + ^“ 3 ^ 

+ Wy [(ct 4 )j^ U + (ag)j^ V + (ag)j^ w] 

+ 0)^ ^ V + (“9)j^ w] 

* {“x * "y "‘'zx'kj - “=xx>jk> 

“z " ^^xy^jk^ 1 ‘ 
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From examining the composition of the inertial force 

note that the first six bracketed terms represent centrifugal 
forces (distance x omega-squared) acting in the deformation co- 
ordinates ^ while the last bracketed terras of Equation 11-97 rep- 
resents Coriolis forces (velocity x omega) . 


MOMENTUM WHEEL COUPLING 


The spacecraft system undergoing analysis may have several ”built 
in” moraentxira wheels* A momentum wheel is generally taken to 
mean a cylindrical or disk-shaped mass that spins about an axis 
that is fixed to a structural hard point of a given body. The 
wheel can be spun up or despun by an electric motor whose rotor 
is part of the rotating mass. The shaft torque that acts to ac- 
celrate the wheel also acts on the body in a negative sense pro- 
viding active attitude control. The shaft torque is generally 
governed by a control law that ^senses" attitude and rate errors 
of the body. In this development a momentum wheel is assumed to 
be inertially symmetric about its spin axis. 



Figure II-4 Typical Body-Momentum Mieel Relationship 



To develop the inertial coupling effects of the typical momentum 
wheel let us consider three unit vector bases: 

[11-98] [e^J “ U, j, kj, 

[11-99] ll^j - [I, m, nj, 

[II-100]and [i j = [I", m', n'J . 
w 

The first triad is the body reference triad for body n, the second 
is a sensor point triad (fixed to point s), and the third triad 
is fixed in the momentum wheel. Now, one of the three unit vec- 
tors of le J is coincident with one of the unit vectors of le j; 

s w 

that is, £, m, or n may be the spin axis depending on the prefer- 
ence of the analyst. In Figure II-4 we have elected to show 

n = n^ as the common, or spin axis. 

The absolute angular velocity of the le J frame can be expressed 

w 

as 

[II-101]w » le J [ R ] ({uj } + {P } 0^ 

w w w s \ s w f 

where is an elementary 3-long position vector (it is null 

except for unity in the first, second, or third locations cor- 
responding to fi., m, or n being the spin axis) and 0 is the rela- 
tive angular speed of the le J frame with respect to the le J 
frame. ^ ® 

With the inertial characteristics assumed (axisymmetry) for the 
wheel, and with the velocity expression of Equation II-lOl the 
total angular momentum vector for the wheel may be written as 

[ 11-102 ]h « le J [J ] {m } 

WWW 

“ [J J ( } + {P } 0 

s w y s ^ / 

with [J^] diagonal with all diagonal values equal to except 

the position corresponding to the spin axis, which is J : J« is 

s X 

the mass moment of inertia about any axis perpendicular to the 
spin axis and is the spin Inertia for the wheel. 
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The torque acting on the wheel (resolved to the [e ] frame) is 

s 

[H-103] T - le^J h 

- N3J ([j„i ti,) + 

- [!>,] [J„l “'v> “) 

where we define an SK* operator such that 


] * SK* {u) or 
s s 


[11-104] 


l— 


— 


— “ 

0 

CJ 

S3 

-UJ ^ 
S2 


^ 1 

Si 

~“s3 

0 

Si 

” SK* 

(U 

S2 

“s2 

S 1 

0 


S3 

— 




— 


The torque acting on body n at point s, due to the wheel is -T 
and it drives the body's quasi-coordinate as 


[ 11 - 105 1 [G'^J - -[OJ 




J (bj {U)„ + (J„l IbJ (U}^ + {P^) J, 6 


-[aj (J.J (».J - la.I [P„l o, 0 


w 


0 


with 


[11-1001 [0,1 - [,\J [^;o|»3] 

and also, as can be easily shown, 

[11-107] [bg] {U}^ - tSK*(j'^R^J [Og] (b^)] [bg] {U}^. 

Wow, the shaft torque is simply the projection of T onto the 
spin axis, or 


[11-103] Tg - [P^J {T} 

■ 0, [1„J [«.! [">„ * 0, [‘■“J ® ■ 

Equations 11-105 and 11-103 allow us to now express the coupled 
equations for body n and several momentum wheels as 
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•„ + 51 + il J„2 iz] 6l P„. 


0 

•^32 


{U} 


n 


J P , bi 
SI wl ‘■ 


J o P o b, 
S2 W2 ^ 


I '^sl 


(G) + [b]"^ {A} 

n n 


u„l [j„l ti,I (U)„ 

0 

+ 

0 

0 


0 


[11-109] + 



(SK*{P }) [b ] {U} J e 
w s ns 


^S2 - ^S2 N t^2l 0n„ 


The inertially coupled body-momentum wheel equations (for two 
wheels) are shown as Equation 11-109 simply for the purpose of 
indicating the form. One may notice that within the equations, 
there effectively resides the original form of the dynamic equi- 
librium equations for body n, namely 


[II-110][m] {U} - {G} + [b]"^ {A} 

n n n n 


which govern in the event that there are no momentum wheels asso- 
ciated with body n. In Equation II-llO we have placed the caret 
(-) over G to represent the right-hand side force vector exclud- 
ing momentum wheel effects. 


Now, on further study of the form of the Equations 11-109, we 
note that if the ^locked” momentum wheel effects are already in- 
cluded in the definition of [ni]^ (which is the standard practice 

when inertially coupling systems together), then the (1, 1) par- 
tition of the coefficients on the left of Equation 11-109 becomes 
simply Also, the second column on the right of Equation 

11-109 is absorbed in having already been accounted for in 

development of dynamic equilibrium equations! 


11-43 



Thus, it follows that in order to implement momentum wheel cou- 
pling with one of the flexible bodies, it is only necessary to 

extend the {U} vector to contain momentum wheel spin values 
n 

(6), to extend the inertia (except for the [1, 1] partition) as 
indicated in Equation 11-109 and to add to the right-hand side 
force vector 

- •'si 'Sll '"n> 

fs, - ^s2 ['■«2j 

The values for shaft torque T that appear in {G } are estab- 

s raw 

lished by a given control law, if the wheels are to be considered 
variable speed* If a given momentum wheel is of constant speed 
(used only for "gyroscopic damping") then the torque equation for 
it is deleted from the form of Equation 11-109; however, its 
effects are still Included in the upper partition of the vector 
{Gn^} (the gyroscopic torque due to constant 0). 

Clearly, the equations of dynamic equilibrium for a body, after 
having been augmented to include momentum wheel coupling, are 
still of the general form 

[II-112]{U}j - + [b]j {A} 


mw 
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F 


GRAVITY GRADIENT EFFECTS 


Attitude dynamics of orbiting spacecraft can be significantly 
influenced by the gravitational force that is distributed accord- 
ing to the system’s position and deformation state. The gravi- 
tational force per unit mass varies (in a central force field) 
simply because different mass particles are at different distances 
from the earth’s mass center. Figure II-5 describes the geometry 
associated with a typical elastic body. 


Earth’s 

Center 



H 

Figure II- 5 Geometry for Gravity Effects on a Typical Body 

For a central force field, the gravitational force per unit mass 
is given as 


[11-113] 


(a 


- GH 
^2 


i 


which, to a first order approximation, is 


Ul-114] ( 

= - g. 



\ 

“/l 

R \ 

c 

1 /j 


where GM is the Earth’s gravitational constant, 
is the typical mass particle, 
g^ is local gravitational acceleration 
is a unit vector directed along 
and c is the origin of the body reference system. 


11-45 



The virtual work due to gravitational force can be written as 



with replaced by differential mass adV. 

The virtual displacement field is expressed in terms of virtual 
displacements of the quasi-coordinates as 

[11“116] Sr ^ Sr +60 x Cpn + "n ) + 6li. 
c c ^ 

In combining Equation 11-115 with Equation. 11-116, the torque 
about point c, due to gravity gradient effects, is 

[II-U7) (tJ . 7^ « s + ^ i-„ « (7 . 

8 c 

where S Is the first mass moment about point c, 

and J is the instantaneous Inertia tensor (deformation dependent) 
for the body. 

The resultant force due to gravity effects is 
® c 

and the force acting in the kth deformation coordinate, is 
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Now, the unit vector has projections onto the body axis system 

that continually vary as the body changes attitude. Let us express 
the unit vector "e^ in terms of direction cosines and the three unit 

vectors associated with the body reference frame as 


[II-120] 




and also define 

[11-121] 


■ =■'* ■ 

[11-122] 

“s - 

W M • 

[11-123] 

w 

= SK* |s| , 

[11-124] 

and 1 

^ fe'l M 


V 


With these definitions and the force and torque expressions of 
Equations 11-117, 11-118, and 11-119, it follows that the first 
three elements of the contribution to the right-hand force vector, 
due to gravity effects are: 


[11-125] 


G 

t gg 


1,2,3 




the second three elements are 


[11-126] 


(g I 

i 4,5,6 


“8c” 





Y i 

I gl 




and the force, due to gravity, acting In the kth deformation mode is 



+ 



(C22) 
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+ (l - j [(b3)^ -b (C33 ),j5j] 

'] 


+ 2y y 

81 S3 


■(bs)^ + (C.3)y £j 


[II-127] * ’'S3 (bs)^ + (C23>^j Ej + (C23>j^ £j 

where the Inertia integrals , (n - l.Z.’-S), and , 

k k j 

(£r,in « 1,2,3), are consistent with the development of Chapter II, 
Section D and Appendix A. 


G. PROVISION FOR INCLUSION OF THERMAL ENVIRONMENTS 


All problems associated with thermally-induced deflections have 
in common the requirement of knowing the spacecraft's attitude 
relative to the sun to deteinnine the effect of solar heating. 

This required information can be extracted, at any point in time, 
from the state vector. It is then necessary to have a model of 
the flexible structure's response, either static or dynamic, to 
solar heating. 

Considerable work has been done on modeling flexible appendages 
in thermal environments* and the results indicate that the response 


*For example, refer to: 

1) Fixler, S. Z-, ’’Effects of Solar Environment and Aerodynamic 
Drag on Structural Booms in Space." J. Spao^araft^ Vol 4, No, 

3, March 1967, 

2) Frisch, H. P., "Coupled Thermally-Induced Transverse Plus 
Torsional Vibrations of a Thin-Walled Cylinder of Open Sections," 
NASA TR R-333, March 1970. 

3) Goldman, R. L., "Influence of Thermal Distortion on the Anomalous 
Behavior of a Gravity Gradient Satellite," AIAA Paper 74-922, 

AIAA Mechanics and Control of Flight Conference, Anaheim, Calif- 
ornia, August 1974, 
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depends on the radiation properties of the booms and the attitude 
relative to the sun. 

The simulation program accounts for time-dependent thermal de- 
formations in the following manner. It is assumed that a model 
exists whereby the structural deformation of a flexible boom (or 
appendage) resulting from solar heating can be determined from 
elements of the state vector and time. This deformation is sub- 
tracted from the actual deformation; the difference is premulti- 
plied by the appendage stiffness matrix. The result is a vector 
of modified, generalized restoring forces for the appendage, which 
is summed into the vector for the appendage body. 

In terms of the development in Sections II. B and II. D where 
-[k]{5} is seen to be the generalized restoring forces (in the 
deformation coordinates) , we note that this is replaced with 
- { 5 ^}). The thermal deformation state } is that 

which must be established from a thermal deformation model. 

In this way, a closed loop response analysis can be achieved 
using external subroutines to develop the thermal deformations . 
Some problems may require only open loop operation if the vari- 
ations of (Cg) in time is slow with respect to general dynamic 

response. 

Rather than building in a rigid (or irrevocable) model of thermal 
deformation, the dynamic simulation program provides the user 
with an interface whereby he can formulate and code a particular 
model, thus latitude with respect to user requirements is retained. 
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III. 


SYNTHESIS AND ANALYSIS OF THE LINEARIZED SYSTEM 


Developments to this point have described the analytical tech- 
niques used to synthesize the nonlinear characteristics of a 
dynamical system consisting of an assembly of interconnected 
flexible (or rigid) bodies. Particular emphasis has been placed 
on spacecraft systems where individual bodies that comprise the 
system may be spinning or nonspinning and may have large excur- 
sions with respect to each other. 

In this chapter we present a comprehensive summary of the tech- 
niques developed for synthesis and analysis of the linearized 
dynamic system with particular emphasis on frequency domain tech- 
niques. For the purposes of this discussion it has been found 
to be convenient to redefine the nature of the system under con- 
sideration and to describe the techniques not in the sense of a 
spacecraft system consisting of interconnected bodies that may 
be subjected to a control system but rather to consider the total 
dynamic system as a pZccnt subject to a oontTolleT, 

Linearization of the nonlinear state equations is necessary to 
permit application of the powerful analytical techniques asso- 
ciated with linear system stability synthesis. Whenever a non- 
linear system can be reduced to a linear system in the vicinity 
of a particular state of interest, it is much more desirable to 
work with the linearized state equations. An additional feature 
related to the linearized system permits the analyst to observe 
linearized perturbation time response characteristics for the 
system. The linearized time response can be quite easily auto- 
mated by recursive formula, which are generally more efficient 
than nonlinear numerical integration algorithms. 



A. 


INTRODUCTORY DISCUSSION 


The mainline nonlinear time domain analysis is structured to as 
senible a collection of interconnected bodies, including a control 
law. The general form of the governing equations may be concisely 
indicated as 


[III-l] = F(Y^,t) i = 1, 2, . . . 

and the form of the function F is the essence of the nonlinear 
time domain solution. In fact, it can be stated that Equation 
III-l is the fundamental basis for the entire DISCOS program. Al- 
gorithms for evaluating the nonlinear state vector time deriva- 
tives (and auxiliary equations) are centered in a subprogram and 
its supporting routines. These same functional algorithms are 
used for linearizing the governing equations about a specified 
state. In addition, it has been found desirable to Introduce 
some new variables Including sensor signals, and control 

torques, B. These new variables extend the nun4>er of equations 
and these additional expressions are linearized along with the 
basic state equations. Additional remarks concerning the use and 
manipulation of the additional variables Is deferred for a later 
section. The remainder of this subsection will address specifics 
relating to the linearization process. 

9 

We first focus our attention on a single variable, y^, and its 

dependence on the system state, Y^, through a known (though pos- 
sibly nonlinear) functional relationship. Arguments begin by 

considering an initial system state, Y^(o), and a functional al- 

* d 

gorithm with which to evaluate the expression, y^^ = ^ 

first express the unknown, y, , in terms of a Taylor s series ex— 

^ i 

pansion about the given state, Y (o) as 


[III-2] y^ = yj^(o) + ^ dY^ + 


a2y. 


i. 


^ dY^dY*" + 


As our Interest lies in the linear part only, the series is trun- 
cated for all partial derivatives greater than one and we have 


[III-3] - yj^(o) 


9Y^ 


dY^ = y 


k.j 


dY 


J- 
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The task at h^d then is to establish the partial derivatives 

indicated as y, , thus yielding an expression of the form (for 
J 

all AY^ « - Y^(o), i = 1, 2, • • •) 

[III-4] AY^ = H AY^ 

1 * J 


Because it would be a nearly impossible (certainly impractical) 

task to generalize determination of the partial derivatives as 

explicit analytical expressions involving the independent state 

variables, we have adopted a numerical approach. This task is 

accomplished by employing numerical perturbation techniques in 

conjunction with quadratic functions to establish the desired 

partial derivatives. Symbolically, we seek to determine the 

elements of H. , such that 

i,j 

[III-5] Y^ - Y^(o) + H AY-^ 

i » J 


where it is assumed that 

• ^ 

1) The functions, Y , are indeed linear sufficiently near the 
state, Y^(o) 

• i 

2) The functions, Y , (although possibly nonlinear) can be rep- 
resented as a quadratic (or lower order) in the neighborhood 

of Y^(o) . 


The basic approach is concisely summarized in two steps: 

1) Establish quadratic coefficients for Y^ in the vicinity of 
the state, Y^(o) 


2 ) 


Evaluate the partial derivatives 

using the quadratic coefficients 
the Independent variables. 


H at the state, Y^(o), 

^ > J 

and perturbation values on 


B. THE LINEARIZATION PROCESS 


With reference to Figure III-l, the quadratic formula can be stated 
in matrix form as 


[HI-6] f(n) - [n2 n ij 
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Figure III-l The Quadratic Formula 

where n is a local special coordinate with origin corresponding 

3f 

to desired to establish the derivative, — , eval- 
uated at • 

In general, the required partial derivative is 


[HI-7] 


3q 


H in 

3n 3q' 


The three values, ^(i+ 2 )^ evaluated via the 

previously discussed functional algorithm, thus these values do 
in fact satisfy Equation III—6. More specifically, consider 
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i+l 

7 

1+2 


\+l 

^ 1+2 


and by matrix manipulation it follows that 


[II1-9] f(n) =• 

""i 

\+i 

1 

1 

-1 

! 

) 

l^i+i ; 


^^^i+2 

\+2 

1 




where the local coordinate, n, is defined to be 
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[III-IO] n 


^ 

^±+2 ~ 


and it can be noted that 

9ri 1 


[III-III n, . O; - i: „ - 


It then follows that 


[III-12] f(n) - n ij 


0 0 1 
1 11 



and if we specify * 1/2 and note that 

we have 


[III-13] f(n) - 



[III-14] f(n) = [n^ n ij 


[III-15] f’(n) 


and, in particular, 

ft I (0) - 


and 


^1+2) ■ ‘’(i) 


^(1) * ^(n=i) 
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1 . 


Comments 


Selection of an initial perturbation value, q(i+2) , from an ini- 
tial specified state, q(o) = » is somewhat arbitrary- A 

value of 1% of the initial value has been successfully used for 
all example problems during the course of the study- In the case 
where the initial value is null, an infinitesimal value must be 
chosen. A value of 1x10“^ has been accommodated in the digital 
code. The intermediate choice of n(i+l) == 1/2 was selected for 
other reasons. Consider first that a single evaluation of a 

3 f 

partial derivative — r is not sufficient to qualify its validity. 

ay 

We have employed an approach whereby two successive evaluations 

of af/ay^ obtained by successively cutting the perturbation in 
half must agree to a predetermined number of significant digits 
(e.g., 5). The choice of n(i+l) = 1/2 requires but a single new 

evaluation for each element in y^ at each successive reduction 
in the perturbation value. In summary, the linearization em- 
ploys an iterative technique to establish the desired partial 
derivatives. 

2. System Resonance Properties 

The linearization process has provided a system of first order 
differential equations that describe the dynamical simulation 
in terms of perturbation variables about an equilibrium state. 

The linearized canonical form appears as 


[III-18] ay’- -H AY^ (i, j = 1, 2, • • ■) 

^ > J 

The coefficients H contain all of the resonance frequency 
^ » J 

properties of the dynamical system. The standard elgensolu— 
tion form is indicated by taking the transform of this expres- 
sion 


[III-19] ^6^-^ 

Extraction of the roots (eigenvalues) from H. then gives the 

^ > j 

roots of the dynamical system. There will be N of these roots 
and any complex roots will appear as conjugate pairs because the 
elements of H. .are all real. The imaginary part of the complex 

^ » j 

pairs represents the resonance (or characteristic) frequencies 
of the system. 
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EXCHANGE OF VARIABLES 


It is often necessary for the analyst to require additional 
variables with which to assess the stability characteristics 
of the dynamical system. These additional variables ordinar- 
ily take the form of plant sensor signals and control system 
output forces and torques. Although the desired variables 
may not be explicitly contained in the system state vector, 

Y , they are known in terms of the state variables through an 
expression of the form 

[III-20] = g(Y^). 

Recall also from previous discussions that either directly or 
through linearization we have established 

[III-21] AY^ =H AY^ 

^ > J 

Now rewriting Equation III-20 In matrix form and Identifying 
variables to retain, Yi, and variables to eliminate, Y 2 , gives 

[III-22] {w[ - [^ 1 1 ^ 2 ] 

and It can readily be established that 
[III-23] |y| = [r] jzf 


where 



and 


1^1 ■{:■}• 

Thus, the state equations for the dynamical system can be 
written (in terms of variables that include the desired plant 
sensor signals and control system forces and torques) as 

Iin-24] iz} = [r]-i |-h^ [r] jzf . 
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and the transformation A = R~^ H R, is commonly referred 

i j ^ > j 

to as a similarity transformation. The matrix is said to 

be the transform of H by the matrix R.* 

1 >3 


The similarity transformation possesses a unique property 

in that the eigenvalues of are equal to the eigenvalues 

of H ! A simple proof establishes this point, 
i >3 

Proof : 


The characteristic matrix of A^j is given by 
[III-25] ^A^j - si ^ = (r~^ j R - sl| =■ R-^ - slj R. 

It follows that Q(s), the characteristic polynomial of , is 

Q(s) » det JA^^ - slj = det R"l|det - sljj det R 

and as (det R~^) = apparent that 

Q(s) = det ^ - slj = P{s) 

where P(s) is the characteristic polynomial of Thus it 

is evident that the matrices and A^j have the same char- 
acteristic equations ’ 


Q(s) - P(s) = 0 

and therefore, the eigenvalues of A, , are equal to the eigen- 
values of H . . 

^>3 

Application of this property now permits isolation of the plant 
and controller, even for a state space representation of an in- 
herently nonlinear system that can be linearized about a spec- 
ified state. Separation of plant and control system variables 
is an important facet of linear system stability synthesis. 


*S. Hovanessian and L. A. Pipes; Digital Computer Methods in 
Engineering, McGraw-Hill Book Company, New York, 1969. 
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1 . 


Evaluation of the Similarity Transformation 



This discussion relates to a procedural approach for determi- 
nation of the similarity transformation matrix, [R], that will 
relieve the user from the burden of having to select those 
variables to eliminate from the original state vector such that 


the auxiliary variables, B and X 


ss 


can become an independent 


constituent of the modified state vector for use in the lin- 
earized studies. With reference to Equation III-22, all of 
the coefficients are known as they have been obtained 

through linearization of the auxiliary equations. The co- 
efficients simply define the dependence of the auxiliary vari- 
ables, w^ , on the original state variables, Y^. In general, 
it is not possible to directly partition the in the Ci and 


C 2 partitions as indicated in Equation III-22, for we have yet 
not made the decision as to which state variables to retain 
and which ones to discard in preference to Introduction of the 


auxiliary variables, w^ . In this light we would like to make 
a best possible choice with regard to which of the variables 

to eliminate from the state vector, Y^, such that the auxili- 
ary variables, w^ , may be included. Many times there will be 
a one to one variable exchange between an element of and an 

element of Y . In any case a variable exchange is necessary 
to structure the total system into the desired plant /controller 
framework whereby the plant and controller can be isolated 
along with the plant sensor signals and the control system in- 
puts • 


The following approach is employed in this simulation to ac- 
complish the desired result; namely, an optimum selection 

from Y as to which variables to eliminate such that w^ can be 
introduced as a part of the state vector. With reference to 
Equation III-22 we can write 


[III-26] 



Our primary focus of attention is now directed to a systematic 
examination of the coefficients such that the variable ex- 
change is accomplished in an optimum manner. We will first 
make note of some size identifications to help clarify the 
discussion . 


III-9 







t-n: 



has size NR by NS 

has size NS by 1 
has size NR by 1 

and 


NJQ = NS + NR. 


Clearly, there exists at least one nonzero element in each row 
of the array. Otherwise Y does not represent an indepen- 
dent set . 

Now a search through the first NS elements of row 1 in the 
matrix array 


[III-27] 



will identify the largest element (absolute value) in row 1. 
Assuming that this element occurs In column JBIG (1 ^ JBIG ^ NS) 
allows us to divide each element of row 1 by this largest 
element and subsequent elementary row operations on rows 1 
through NR will eliminate those elements below the pivotal 
element in column JBIG* 


This procedure is repeated for each of the NR rows contained 
in the matrix and the following observations are noted; 

1) the appearance of a one (1.0) in a row identifies a vari- 
able that will be eliminated in preference to inclusion 

of an element of w^ . 


2) The absence of a zero or one in columns of a given row 
indicates which variables will survive the exchange process. 

3) All variables in (NR of them) will become part of a new 
and independent state vector (the modified state vector). 

4) The transformation, = 1 • • * NS) can be con- 

structed from the matrix that remains after the procedural 
approach has exhausted all of the NR rows of the expression 
III-27. 


2. Illustrative Example 

A simple example is presented to further describe the actual 
mechanics used for evaluation of the similarity transformation. 
Linearization of the auxiliary equations has established 
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and to implement the proposed algorithm, we first write down 
the C matrix augmented with and further identify the elements 
of C as 


I til ^fbiibi2bi3bi4 -1 "I 

I Ij’2lb22b23b24 “ij 


For Illustrative purposes assume that |bi 3 l > l(bii, bi 2 » ^ 14 )! 
and hence bi 3 is the pivotal element for the first row. It 
follows then that row 2 is modified by the relation 


ele2j = - b23 * bx j + b2j 

t>is 


giving the matrix 


bll 

bl2 ^ 

^13 

1 

t>14 / 
'^13 

1 

bi3 

0 

-b2 3 

bl 3 

_ +^2 1 

~t>2 3 

bl3 

+t>22 

0 

-W3 

+b24 

b23 

bis 

-1 


Cll 

C 12 

1 

Ci4 

CIS 

0~ 

^21 

C 22 

0 

C24 

C25 

-1 


We continue the process by first dividing row 2 by the largest 
element and again using row operations to eliminate the other 
element in that column. Again, for illustrative purposes, we 
will assume we know the pivotal element of row 2 (say C 22 ) 
row 1 will be modified as follows by row operations 


ele^j = 
giving 


-Cl 2 




'22 


-C12 

4.^ ^^22 
+Cl 1 

0 

1 

-C12 ^ 
. C 22 

+C14 

C2 S 

-Cl 2 

^22 

+C15 

C12 

C 2 I 



C24 

C25 

1 


1 

0 





C 22 



C 22 

C 22 

C 22 


dn 

0 

1 

di4 

di5 


d23 

1 

0 

d24 

d 2 5 

d26 


III-ll 



from which we establish the desired similarity transformation as 


” 1 

0 

0 

0 

-dai 

-dji^ 

-d25 

-d26 

-dll 

-dik 

-di5 

-di6 

0 

1 

0 

0 




Now it follows that the original state variables are written in 
terms of the modified state variables as 


A 

1 

0 

0 

0 ' 

pl\ 

IY2/ 

-d2i 

-dan 

-d2 5 

-d26 

\yi, I 

1 M ' 

-dll 


-dx5 

-di6 

1 M 

1 Y4I 

_ 0 

1 

0 

0 

( ^^2 1 
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SYSTEM TRANSFER FUNCTIONS 


The entire system transfer function synthesis can be concisely 
summarized in a chronological sequence of steps that began with 
linearization of the coupled me chan leal/ control law equations 
that govern the dynamical motion; This process included linear- 
ization of additional equations that contained specific variables 
required for further consideration in the' stability analysis; 
namely, plant sensor signals and control system outputs • A 
similarity transformation has been introduced which in effect, 
exchanges original state variables for these desired sensor sig- 
nals and controller outputs such that the resulting modified 
state vector still is representative of an independent set of 
state variables. The resulting system of state space equations 
is later identified as Equation III-28. 


The system characteristic matrix, j > provides the basis for 

evaluating the coupled mechanical/control system resonant char- 
acteristics (natural frequencies) as well as providing the fun- 
damental basis for specification and determination of the various 
types of transfer functions. The next subsection addresses some 
of the more specific details regarding specific transfer function 
relationships- A particular transfer function is identified by 
a type along with the desired output /input variable designa- 
tion. An eigenvalue problem is then stated, which leads to 
determination of the numerator roots (zeros) and denominator 
roots (poles) for the particular transfer function. Once the 
poles and zeros are known for a transfer function, this infor- 
mation can be further processed and displayed by any of the 
conventional display modes; Bode, Nichols, Nyquist, and/or 
root locus- 



The conventional block diagram representation for the coupled 
plant/controller system (Figure III-2) provides additional in- 
sight for determination of system transfer functions. 



Figure Flant/Contr oiler Block Diagram 

The first-order differential equations for the system are written 
as 

[III-28] = A , + Bg 

ij ij 

and it is helpful at this point to express the equation in m- 

trix form and indicate the separate partitioned subsets of Z , 

A . B , R„^ , B and R ^ as 

*ij’ ’ '^ij ^ ®ij ® 



The following observations are noted; 

^31 “ ^ '^T^ “ *^s^ ° ^ 

aiti = 0 b^2 “ “^24 bg2 = Q 

ai3 = 0 b^3 “ 0 b^3 = 332 

a23 “ 0 b^4 - 0 b^4 - 342 
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and Equation III-29 can be restated as 



ail 

ai2 


ai4 

a2i 

322 


324 


332 

^33 

^34 


342 

343 

a44_ 




Equations III-30 are the operating basis for stating particular 
transfer function relationships for the plant/controller system. 


The general procedure is to establish a system transfer function 

between inputs iL, and R and outputs X and B. Loops may be 
Is s s 

opened to provide open loop information by manipulation of the 
coefficients to prohibit certain feedbacks. 

To symbolically describe specification of a transfer function 
we begin by consolidating the b coefficients and taking ttie La^- 
place transform of Equation III-30 to give 

[in-31] [isj - [a] + [b] 


or 


[III-32] 


['“] - Hi - H {“(s)} 


and then employ Cramer’s Rule to evaluate a given element 


due to a particular input U 


[III-33] Z 


, X /U, . 

(s) (s) 


aug Is - A 


(s) 


where 


Is - A 


and where aug [is - a| is accomplished by placing column q of b 
into column p of jis - A ]• 

The Q— R algorithm is a useful tool with which to extract the 
indicated determinants in Equation III-33* 


J. G. F. Francis, "The QR-Transformatlon - A Unitary Analogue 
to the LR-Transformation." The Computer Journal, Volume 4, Octo- 
ber 1961 (Part 1) and Volume 5, January 1962 OPart 2). 
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1 . 


The Root Extraction Process 


With reference to Equation III-33 it is desired to evaluate both 
the numeratoxf and denominator roots. The denominator root extrac- 
tion is straightforward in that we wish to find pi, P2> P3» *** 
from an expression of the form 

D(s) * det^[I]s ~ [A]) 

such that 

[III- 34] D(s) - (s - pi) (s - P 2 > ••* (s - P^) 

This evaluation is completed by extracting the characteristic 


n (s - pj. 

i=l 


roots of the matrix . In general these roots will be complex 
because is not symmetric. 

The process employed for evaluating the numerator is best illus- 
trated with an example. Consider that we have the (4x4) charac- 
teristic system matrix, 

[Ayl 


'ail 

ai2 

ai3 

314 

^21 

322 

^23 

324 

aai 

332 

^33 

^34 

am 

342 

343 

344 


and the column of coefficients b^ which premultiply the desired 

input variable Further, let it be desired to obtain the 

transfer function relating output of the third variable in the 

state equations ys to the input U^. 

The state equations for this system would appear as 


[III-35] 



■311 

^12 

^13 

314" 

321 

322 

323 

324 

331 

332 

333 

334 

.341 

342 

343 

344. 



and, with reference to Equation III-33, the numerator is 
N(s) = aug|ls - a| or 
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After performing elementary row operations. Equation III-36 can 
be restated in the form 


s -an + a3i bi/bs 

-aij + a32 bi/b3 

-ai4 + 334 bi/b3 

[m~37] N(s) - b3det -azi +^31 b2/b3 

8-a22+asi b2/b3 

-324 +334 t>2/b3 

-ai,i +331 l>4/b3 

-342 +^32 t>4/b3 

8-344+aS4b4/b3 



or, in symbolic terms as 
[HI-33] H(s) - b 3 det|[Is] - [a] [ 

where the matrix a is given as 


an -331 bi/b3 

ai2 -332 ti/b3 

314 “®34 bi/b3 

^21 -331 b2/b3 

^22 ~a32 ^ 2/^3 

324 -334 b2/b3 

341 “331 t»4/b3 

342 -332 t>4/b3 

a44 -a34 b4/b3 


Note that the previous expression for N (s| is finite only if b 3 ^ 0 
and the question is — can b 3 realistically be null? The answer is 
yes as the following example indicates. 

Example 

Consider the simple mechanical system consisting of two masses 
connected by a single spring/dashpot combination as shown in the 
sketch. 



The state space representation is 
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mi-- 



-c/mi 

c/mi 

-k/mj 

k/mi 


Vnii 0 

1 ^2\ “ 

cimz 

-c/mz 

k/m2 

-k/m2 

j X 2 I + 

0 Vmz 

jxi| 

1 




1 M 

0 0 



1 



(xz) 

0 0 


1^1 

1^2 


and the frequency domain (or transformed) equations In s are 
[[I]8 - [A]]|^Xi(s)'j - p/mi 0 IjFij 


/"Xi(s)\ = 

' /m^ 0 

1 X 2 ( 3 )/ 

0 ^/vnz 

JXi(s)j 

0 0 

(x2(s)) 

_0 0 


where 

[A] - r “c/mi c/mi -k/mj k/mi 

c/m2 -c/m2 k/m2 — k/m2 

10 0 0 

0 1 0 0 _ 

Consider now the transfer function Xi (s)/Fi where the augmented 
numerator is 

det I Vmj^ -c/mi k/m^ -k/m^ 

0 8 + c/xa.2 -k/ma kMj 

0 0 s 0 

0-1 0 8 

and the pivot element is the (1, 1) element or Vmi fO. On the 
other hand, the transfer function Xi(s)/Fi has the augmented 
numerator 


“(») 


N(a) 


det 


s + c/mi 

-c/mi 

Vmi 

-k/mi 

-c/ia2 

s + c/m 2 

0 

k/m2 

-1 

0 

0 

0 

0 

-1 

0 

s 

3t element 

is the (3 

, 3) 

element » 
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The problem we now address involves evaluation of the nximerator 
determinant N(s) when the pivotal element is null. The particu- 
lar mathematical problem may be restated as 

[III-39] n(s) - det I [i]s - 
where 




is the identity matrix 
element. 



of size N with a null diagonal 


[i]x (w, 


Addition and subtraction of the quantity 
bitrary constant not equal to one of the roots of 

[III-40] N(s) - det 

J ^ - 

and if we define (s - x) = Vp> there results 

[III-41] 


where x is an ar- 


[i] (s - X) - [[i] - [i] x] I 

line (s - x) = Vp> there resu 
N(s) =■ (- 1 )^ det I A - ix j|det j Ip - (a - lx) I ^ 


yields 


The roots, (p^, i=l,Nl are found as the eigenvalues of 
pression 

(UI-42) [[a] - [i] x] [l] 

and the eigensolutioh permits n(s) to be written as 

[III-43] N(s) = ("iV^ det I A - lx 
N 

P I 



the ex- 



We now make the following observation: a p^ equal to zero im- 
plies a root at infinity (or a characteristic polynomlnal having 
order less than N) . Thus, the null p^’s are eliminated from the 

expression giving the characteristic poljnninal an order n, which 
is less than N. It is a rather common occurence for the number 

of zeros (order of N(s) ) to be significantly less than the num- 
ber of poles (order of D(s) )• With reference to Equation III- 
43, the numerator expression, N(s) can be written as 
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[III-44] U(s) - (-1)^ det 


A - ix| 

(d - Pi ) (1 - ) ••• (1 - Pn 

■ 

t " " M 


and, recalling that p =* > yields 

S“X 


[III-45] n(s) » (-1) det A - Ix 
n / V I 

n 

i“l 

Next, we note that 

n ^ / -1 \ 

[ni-46] n (X - s ) “ /7 77 

i-1 ^ i=l\ ! 


(s - Sj) (s - S 2 ) ••• (s - s^) 


and it follows that 

[III-47] N(s) - (-1)^ /7 Pj 

i=l 


A - Ix 


n 


/7 (s - s^), 
1=1 


The numerator root gain, 
[III-48] - (-1)'^“" ^1 


k-, , can now be identified as 

A-ixl 


and the Bode gain, k^, for the ninnerator is 
m 

k_ * k_ (-1)^/7 s where m ^ n. 

B ^ '4 1 

2. Transfer Function Classification 


With reference to Figure III“2 it is possible to directly iden- 
tify six transfer function types* Each type is characterized by 
the specific variables involved and by the presence of feedback. 
Additionally, a seventh type will also be described whereby cer- 
tain of the control variables feed back and others do not. This 
type is similar to an open loop transfer function but treats se- 
lected channels of the controller as part of the mechanical sys- 
tem (plant) . During the course of this discussion it will be- 
come apparent that additional transfer function types are easily 
accommodated by rather simple raanipultaions with the system char- 


acteristic matrix. 
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In general it should be noted that the process of obtaining the 

desired transfer function involves but a few basic steps. The 

transfer function characteristic matrix, fR. and the desired 

ij 

force coefficient vector, are obtained directly from the sys- 
tem characteristic matrix A 


ij* 


These two matrices are then put 


in a form such that the Q— R algorithm can be employed to extract 
system roots. 

Type I (Plant Only) 

Type I is the forward path transfer function for the plant with 
no feedback and is of the form 


[HI-49] X^P/R^ - G(s) 


The control variables 6 and control outputs, B , do not feed 
back into the plant. The matrix expression depicting the system 
of Interest is 



ail ^12 

a2i a22 



The matrix, , to use in the general expression given as Equa- 

tion III-33 is referred to as |R .. or the reduced A . . matrix, 

iJ ij 


[in- 51 ] IR 


ij 


an ai2 
321 ^22 


The augmented matrix is obtained by removing the colxamn cor- 
responding to the input variable, from the expression b^ and 
inserting this column Into the column in |Rj|^ j » which corresponds 
to the desired output, X The resulting transfer function is 

3 S 

then given as 


[HI-521 X P/d - J j 

^ 1 1= - IR| 
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Type II (Controller Only) 

Type II represents the feedback path, H(s) , for the controller 
only. The desired transfer function relates control system out- 
puts to sensor signal inputs, 


[HI-53] - h(s). 

The reduced characteristic matrix corresponding in- 

put coefficients, given as 


[III-54] 

"*33 *34 

*^k = 

1 

CM 

m 

*43 *44 

— 

*42j 

Type III 

(Open Loop, GH) 



Type III falls within the framework of the classical open-loop 
transfer function designation and relates control system out- 
puts to external plant inputs K^. The algebraic expression for 

a given output variable, due to an external input, R^, is in- 

dicated as 


[III-55] B^/R^ - (GH) 

The open-loop system characteristic matrix 
input coef f iciencts, b^^, are 



and corresponding 


[III-56] 


r*ii 

*12 



*21 

*22 




*32 

*33 

a4 3 


*42 

*43 

*44 


b,, 

ik 


-am . 

"*24 

0 

0 


Previously it was noted that 331 = a^ » ai 3 “ 823 “ 0 
dition, the partitions am and 324 to zero to prohibit 

the B^ feedback. Thus, the loop is opened to establish GH, the 
open-loop transfer function in s. Note that the negative sign 

in the b^, coefficients simply indicates that the B feedback is 

Tjj 

negative with respect to the external plant inputs, R^. 
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Type IV (Open Loop> HG) 


An additional open-loop transfer function is often desired to 
assess the plant sensor signal outputs due to controller noise in- 
puts. The transfer function then relates sensor signal outputs, 


to control system noise inputs, 


The plant sensor sig- 
nal vector does not feed back into the system so that we have 


[III-57] X - (hg) , - 

ss s ' ' (s) 

and the system characteristic matrix, IR 
put coefficients, b^^, are identified as 


ij* 


and the external in- 


% 


[III-53] 


ail 

ai2 


ai4 

^21 

322 


324 



^33 

^34 



^43 

344 


'ik 


0 

0 

^32 

^42 


Note that the a^z ^42 partitions have been nulled to elimi- 
nate sensor signal feedback. 

Type V (Closed Loop - Control Ratio) 

The system control ratio is given as the transfer function that 
relates plant variable outputs to externally applied plant inputs 
with the control system entirely active. We express this transfer 
function as 


[III-59] X P/r2 

s s 


\l+GH/(s) 


and the system characteristic matrix and the external input 


coefficients b,, are identified as 
ik 


IR 


ij 


[III-60] 


”ail 

ai2 


314 

321 

322 


324 


332 

333 

334 


342 

343 

344 

m 


ik 


-ai4 

-324 

0 

0 


The negative sign in the matrix indicates that the feedback is 
negative. 
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Type VI (Closed Loop) 


An additional closed-loop transfer function has been accommodated 
within the digital simulation. Specifically, Type VI relates 
plant sensor signal outputs to sensor signal noise inputs with 
all control system loops active. The transfer function is sym- 
bolically indicated as: 

[III-61] X ^ “ (transfer function) 

s s s 

where the system characteristic matrix, corresponding 

input coefficients are identified as 

IR 


'id 


[III-62] 


an ai2 


314 

• ^ik ■ 

“0 

321 ^22 


324 


0 

332 

333 

^34 


as2 

342 

aii3 

344 


342 


Type VII (Quasi-Open Loop) 

An additional transfer function type is identified here and re- 
ferred to as quasi-open loop. It is of the open loop type in 

that we are interested in control system outputs, due to plant 

variable inputs, R^. For example, suppose that for a multi-channel 

control system (such as azimuth and elevation), we desire outputs 

on the controller channel that do not feed back and that the 
other channel is active in that it feeds back into the plant. 



For the configuration indicated, a typical Type VII transfer func- 
tion (TF) would be given by 

B 2 “ (transfer function) 
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and the form of the system characteristic matrix, 
input coefficient matrix would be 


IRjLj 

[III-63] 





— 

an 

ai 2 


314 

321 

a 22 


324 


as 2 

a33 

334 


342 

a43 

344 


'ik 


-ai4 . 

-a24 

0 

0 



and plant 


The subpartitions a^ and a 24 indicate modification of the origi- 
nal partitions ai 4 and a 24 « Specifically, a^ is a subset of a^j 

obtained by keeping only those n columns of a^ that correspond 

to the variables that feed back to the plant* 


3* Transfer Functions - Polynominal Description 

This subsection is addressed to implementation of control system 
transfer functions described as the ratio of two polynominals in 
the frequency domain, s* Specifically, we consider 

[III-64] TF - P(s) /q(s) 

where 

Q(s) » ao + ajs + 328^ + 338 ^ + ••• + 

and 

P|s) » bQ + b^s + b23^ 4“ * • • + 

Because the previously described governing equations have been 
stated in cononical first-order form, we propose to restate the 
polynominal description for the transfer function in the form 

[HI-65] 6^ + B^U 

The block diagram for the system is 


u 


6^ 


Q(«) 



[III-66] from which we write 



U 
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and expansion of the implied operator in s results in a differ- 
ential equation of the form 

n n-1 , m m— 1 

[III-67] a 6 + a - 6 + ••• + ai 6 + an6 «b U + b - U + ••• + bi U + bnU 

^ n n-1 ^ ^ m m-1 ^ ^ 

wnere 6 - d 6 , 



In general, the order of P |s) will be no greater than the order of 
Q (s) or m^n. 

a, m«n 

We divide Equation III -67 by a^ to obtain 

n n-1 • ni m-1 . 

[III- 63 ] 6 + 6 + Cq<S = d^ U + U + ••• + diU + doU 

b^ 


where C. 


— and d , 
a i 

n 


An example will be used for illustration* 

Example: Consider the equation with m=n= 4 , 

♦••• ••• •• * ■#*• *•• • 

6 + C3 6 + C2 6 + C16 + Cq(S - di4**U + d3*u’ + daU* + djU + dflU 


or. In operator form 

s‘*iS + a^ C36 + C26 + Cis6 + Co<S •• s‘*d4U + s^d3U + s^d2U + s diU + doU. 

This can be rewritten as 

s‘*(6 - di^u) + s3 (c3& - d3U j + s 2^C2<5 - d2U ) + s |Ci 6 - diU j + (^Cq6 - doU)= 0 

and the substitution 
6i = 6 - di+U 

permits a reduction in order to 

+ C36 - d3uj + “ ^ 2 ^) s|Ci6 - d^uj + 

we can again introduce a new variable 
^2 * (^1 ^ 3 ^ “ ^3^ ) 
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and rewrite the previous as 

(^2 8|ci5 — diuj + |Cq6 “ <^0^ ) “* 

It follows that if we define 

m 

^3 “ ^2 ^2^ “ d2U 

there results 

3(63 + C16 - dju) + Co<S - dgU * 0 , 

and the substitution 
64 * 63 4 * “ diU 


gives 

64 “ “Cq 6 + doU» 

The variable 6 can now be eliminated from each of the above ex- 

'th 

pressions and the results generalized to n order systems. 


The result is concisely stated as a matrix equation that is recog- 
nized to be of the desired form initially given as Equation III-65, 

U 


[III-69] ( • 




— 





— 





"\ 




1 

0 

• 

• 

0 





- C -d 


n-1 








n-1 

n-1 n 

h 



0 

1 

• 

• 

0 


^2 


^ n 

- C .d 



n-2 









n-2 

n-2 n 


a 

• 

m 

• 



• 


• 

“h 


« 



• 

• 

m 



• 

< 

« 



• 

• 

6 - 


-Cl 

0 

0 



1 


6 , 


di - 

C 1 d 

n-1 








n-1 


^ n 

5 

n 


-Co 

0 

0 

» 

a 

0 


6 

n 


do - 

Cod„ 



- 





— 





J 


where and 5, the original variable of the equation, are rela- 
ted as shown previously and U is the input variable to the trans- 
fer function expression as Indicated in Equation III-65, 


b, m < n 


The general expression for the case where m<n is easily accomo- 
dated by restricting the coefficients to reflect the limit m. 

Commonly, only the dQ coefficient will be finite. 
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4, Frequency Response 

Transfer function poles , zeros ^ and root gein can be converted to 
the standard Bode form for frequency response by combining time 
constants, damping, and resonant frequencies as 

N2 


[III-70] TF » kg 




Ml / \ 

\ “J "3 7 


where the Bode gain Is 


n 

n\ 

i=l 

m 


where k * root gain and 


X = system constants 
? = system damping at frequency u) 
w = system resonant frequency. 

T]xe frequency response is then calculated by substituting Jw 
s and evaluating the transfer function expression at various oj's. 
The digital simulation uses a vernier frequency incrementing ap- 
proach that automatically Introduces smaller frequency increments 
near the poles and zeros. This variable frequency Incrementing 
technique permits better transfer function resolution near the 
resonances where amplitude and phase can vary rapidly. 


5. Root Locus 

The root locus method of analysis and design is based on the re- 
lationship between the poles and zeros of the closed loop transfer 
function and those of the open loop transfer function. The method 
is used to determine the location of the roots of the character- 
istic equation as a function of a single open loop gain param- 
eter, The locations of these roots are indicative of the relative 
system stability. The analyst may use the method as a design tool 
by adjusting the poles and zeros and the open-loop gain parameters 
in such a way as to yield a closed loop system with satisfactory 
critical frequencies (poles and zeros). 
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To further describe the theoretical basis for the method we re- 
fer to the conventional control ratio for a feedback system as 
shown in Figure III-3, 



Figuj^e III~Z Conventional Feedback Control System 

The control ratio C|s)/R(s) is 

[III-71] _ g(s) 

R(s) 1 + G(s) h(s) 

and the open loop transfer function is identified as 

a ratio of two functions in s, 

[III-72) G(s} H(s) - k ’ 

Q(3) 


The characteristic system equation is 


[III-73] 1 + G(s) h(s) = 


or 

1 + k 


Zisi 

Q(s) 


0 . 


0 


The conventional root locus plot portrays the loci of the values 
of s that satisfy the characteristic equation as k varies from 
zero to infinity and we note 

1) at k=»0, the roots of the characteristic equation are equal to 
the roots of Qfs), which are the same as the poles of the open 
loop transfer function, k p(s ) ; 

Q(s) 
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2) as k approaches infinity, the roots approach the roots of P(s), 
the open loop zeros. 

Thus, as k varies from 0 to infinity, the loci of the closed loop 
poles migrate from the open loop poles to the open loop zeros 
and the direction of migration depends on the sign of the open 
loop gain parameter, k. 


Rewritting Equation III-73 yields a more conventional expression 
for the characteristic equation as 


-1 


[HI-74] k P(s) 

Q(s) 

and two conditions are required; 


1 ) 


k P (s ) 


1 ; 


2 ) 


’Q(3) 

/v{a] /q(s) = 


130°, k>0 


The first of these conditions can be expressed as 
l Q(s) | 

P(s) 


for those values of s that satisfy the angle criterion. The con- 
ditions that govern the migration of the roots in the complex 
plane can be solved by an iterative procedure. The iterative 
procedure for evaluation of a single root locus* is described in 
Appendix E. 


E. LINEAR TIME DOMAIN RESPONSE 

The linearized canonical first-order system of equations can also 
provide a basis for studying system time history in terms of per- 
turbations about a specified state when the system indeed behaves 


*Welch, Raymond V. National Aeronautics and Space Administration 
Goddard Space Flight Center, Branch Report No. 254, October 2, 
1973. 
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in a linear manner in the vicinity of the state. The nonhomogen- 
eous form of the equations was the basis for determination of 
system transfer functions and appeared previously as 

UlI-75] Z^ + u’*'(t). 

The external system inputs are the elements of U . It is con- 
venient to establish the solution for the above system through 
use of a recursive formula numerical integration procedure rather 
than through the Runge-Kutta approach. 

Consider the Adams’ corrector formula* at time t+1, 

[ln-76] ■ "t * p \+l + 19 - 5 

where h is the incremented time step. 


Application of this formula to our system of equations gives 


=» 2^ + 9 A 

' t+l t 24 ^ t 


t+l + 9 + 19 - 5 


and manipulation yields the solution for all the z at time step, 
t+l 


izi 

II t+l 




I 

(t+l 


+ 19 { Zf . - 5 |z 


•lt-l 




Note the requirement for z at time step t-2; hence, the require- 
ment for a starter (e.g., Runge-Kutta) to initiate the solution 
process. 


F. Schied. "Theory and Problems of Numerical Analysis." Sahawn'a 
Outline Series, McGraw-Hill Book Company, New York 1968. 
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APPENDIX A-^DEVELOPMENT OF THE INERTIAL INTEGRALS 


In the development of the equations of motion (Refer to Chapter 
II, Sections B and D.) there are certain inertial integrals 
identified that are required to account for the deformation-depen- 
dent inertia matrix and that are Involved in calculating the ef- 
fects of centrifugal and Coriolis forces. 


[A-1] 


The basis for calculating these integrals is a triple matrix pro- 
duct involving a so-called discrete mass matrix [M], which is as- 
sembled by use of finite element techniques, and which may be used 
in calculation of vibration modes. The other constituent of the 
triple matrix product is a modal transformation that transforms 
ordinary velocities, associated with the finite element model, to 
the velocities of the vector. 


Let us refer to the transformation as [<(>], thus the triple matrix 
product is 


[m] - 


which is the basis of the kinetic energy expression of Equation 
11-21. Now, the mass matrix [M] is invariable with respect to 
the body’s deformation. The modal transformation [4>] does, how- 
ever, depend on the {?} in a linear fashion, or we may expand 
l^] as 

[A-2] [<()] ■ 4- [A4>] , 

with [4“]^ a matrix of constant elements and [A((y] variable with 
respect to deformation. 


On substituting Equation A-2 into Equation A-1 and referring to 
Equation 11-86 it follows that 


[A-3] [mj - 

[A-4] Cj = 


and 




A-1 


[A-5] 



[A-6] 


[A-7] 


Assume that the finite element model of the body has a ^’global” 
cartesian frame in which the ordinary velocities are measured, 
and further assume that the generalized coordinates of the finite 
element model are grouped (or ordered) such that all the x-trans- 
lations are together, followed by all the y- then z- translations, 
and that the translations are followed by sets of x, y, and z ro- 
tations, With this implied ordering, it follows that the discrete 
mass matrix is partitioned in the form: 


m 

XX 

m 

xy 

m 

xz 

m 

xp 

m 

xq 

m 

xr 


m 

yy 

m 

yz 

m 

yp 

m 

yq 

m 

yr 



m 

zz 

m 

zp 

m . 
zq 

m 

zr 




m 

PP 

m 

pq 

m 

pr 

(S 

YMMET 

RIC) 


m 

qq 

m 

q^^ 






m 

rr 


with p, q, and r corresponding to rotation coordinates about x, 
y, and z axes, respectively. Similarly, the modal transformation 
is partitioned as 





{z+n^} 

-{y+riy} 

{1} 




-{z+nz> 




{1} 


f — 1 






{1} 


{1} 








{1} 








{1} 






Each square subpartition of Equation A-6 has rows equal to the 
number of structural joints (collocation points) of the finite 
element model, as does each subpartition of Equation A-7, The 
submatrices in the last column partition of Equation A-7, 

([h ], [h ], [cJ ])i have columns equal to the number of 

deformation modes used to represent the body and are matrices 
of modal translation and rotation amplitudes. 
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[A-8] 


[A-9] 


[A-10] 


[A-11] 


The form of and of [A^i] Is seen immediately from Equation 

A-7 in that the only nonzero parts of lA<J)] are due to the {n) 
vectors. The [(j)] matrix is effectively a kinematic velocity trans- 
formation consistent with the form of Equation H-25, and it 
follows that 


V ■ tV"*’ 

<"y) - IfylU), 

and {q } - [h ]{^} • 
z z 

In the Equation A-4, there Is seen the product of two constant 
matrices*, namely The two triple products on the right 

of Equation A-4 require evaluation of only the first three row 

partitions of [M] [(^ ] . Thus let us define 

o 

(The first 3 row partitions of [M][<()^]) * 


1^x1 1 

1^*2 1 


|^4j 

l^s! 

^xej 

' -V ' 
1 


i 


1^3} 

fv! 

j^sl 

jv| 

V 


_1M i 



I^zaI 

!^z 5 j 


*"zk' 



with, for example, 

■ [”xp] 'y’ - ["xyj 

and 

It is unnecessary to expand each partition of Equation A-9; the 

product is numerically obtained and the examples of Equa- 
tions A-10 and A-11 are just for purposes of illustration. 
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Now with reference to the intermediate constant matrices given 
by Equation A-9 and the definitions of Equations A-4 and A-5, the 
following Inertial integrals are developed (the reader is urged 
to refer back to Chapter II, Section D, particularly Equations 
11-88 and 11-89): 

|“>J ■ f-y] - [^4 [»z 

[A-13] . jp^jJ [hj - 

[A-14] a, =■ P. M - P.[^' 

^ ZD y yo 2 

[A-15] 04 - 

[A-161 |,5j - |^P^ 3 j 

[A-17] |«sj - |h,] - i^p^sj k' 

[A-18] aj =■ P, [hi - P, fi* 

' y4 X x4 y 

[A-19] ao * P c M “ [P c: 

° y5 X x5 y 

[A-20] 

‘ ^ y6 X xo y 

[A-211 |b,j - [p^jj - [py,j [h^- 

iA -221 - jp^jj - |p^2j k 

[A-23] L - Py3 y - p^jik 

[A-241 |b4j - - [P^iJ [hz| + I^PyjJ - I^P^jJ jhy 

[A-25] [bsj - [P^iJ - jVyiJ k] + ['■y3j k’ - [^ 3 ] k 

[A-261 jbej - [P^jJ k f*:] + [^ 3 ] k] - k3j fz 

N' klk-'l ■ 

n-nm-nn 

(A-30) jb,,]. - klkylf"^ 


A-4 





> 

I 

Ln 



APPENDIX B--DEVELOPMENT OF ROTATION TRANSFORMATIONS 



There are 12 possible orthonormal rotation transformations, in 
terms of Euler angles, that the analyst may choose from in order 
to orient one orthogonal triad with respect to another. For each 
one of the 12 orthonormal rotation transformations there is an 
associated rotation transformation that is not orthonormal and 
that is used to transform angular velocity projections (onto a 
nonorthogonal vector basis) , which are time derivatives of Euler 
angles, to projections (onto an orthogonal vector basis) that 
are commonly referred to as time derivatives of angular quasi^ 

coordinates (a> , w and o) ) . 

X y z 

It is possible, for purposes of digital computation, to automate 
the generation of these transformations, given a selected order 
of rotation. It is the purpose of this appendix to indicate the 
steps and numerical manipulations that are required,. To this end, 
let us consider one of the 12 types (say a 2 - 3-2 permutation) as 
an illustrative example. 


Consider the two orthogonal vector bases, whose relative orienta- 
tion we want to describe, to be 


[B-1] {a} 


and 


I 

J 

K 


[B-2] 


{e} 


Now if ©i, 02 , and 63 are the three successive Euler rotations 
about axes (2-3-2) respectively, then it follows that 


{a} - 





[B-3] 

COS0 1 


sin 0 1 

i*' 



1 


j " 


-sin9i 


COS01 

k" 


B-1 



[B- 4 ] 


{e'} = [T2]{e"} 


= 

COS02 

-sln02 


1" 


sin02 

cos 02 


J" 


. 


1 

ic" 


and 

{I"} - [T3]{I} 


[B- 5 ] 


COS03 

1 

sin03 

■ ~ 
i 

I 



-sln03 


COS0 3 



On combining Equations B- 3 , B- 4 , and B -5 there results 

[B-6] (I) - iTiJETaJtTgKl}. 

Now, a 2 - 3-2 permutation meOTs that the first rotation (Sj) is 
about the 2nd axis of the {a_} basis, the second rotation (62) is 
about the 3rd axis of the {e^^basis and the third rotation (63) 
is about the 2nd axis of the {e**} basis . 

Consider the following reference table, which shows the correla- 
tion between Euler rotations and the corresponding axis: 


Table B-1 Correlation of Euler Rotations and Axes 


Type 

m 

2 

3 

m 

5 

6 

Bl 

8 

9 

10 

11 

12 

0j about 

n 

Bl 

Bl 

Ha 


ESI 



3,K 

EEl 

3.K 

3,K 

02 about 

2J' 

BQ 

ijj 

Vl 


3,k' 


m 

la 

la 

D 

2,r 

2,r 

63 about ' 

HBl 

Ha 

la 

CM 

m 

HBl 

2,j" 

m 

2,j" 

3,k" 

Bl 

— tt 

l,i 


Now it is clear that the elementary rotation transformations ([Tj], 
[T2], and [T3]) always Involve 0 i, 62 » respectively, but any one 
of them may have three different forms depending on the axis asso- 
ciated with its rotation. That is, when 0 ^ (i ■ 1 , 2 , 3 ) is about 

axis (1) then 


[B- 7 ] 

[T^] - 

'1 

COS0^ 

-sln0^ 



sin0^ 

COS0^ 


when 0^ is about axis (2) , 


B -2 































[B-8] 


I 


[B-9] 


[B-10] 

[B-Il] 


[B-12] 


[B-13] 


[T^] 


- 

cos9^ 

1 

sin0^ 


-sinOj^ 


COS0^ 




- 


and finally, when 0^ is about axis (3), 


[T^] 


COS0^ 

-sin0^ 


sln0^ 

COS0^ 

1 


Thus, it is evident that one need only specify a rotation type 
(referring to Table B-1) and the three Euler rotations to create 
his required orthonormal rotation transformation Equation B-6. 

The associated rotational velocity transformations are developed 
as follows. Consider, again, the 2-3-2 permutation. For thi^s 
case, it is possible to express the angular velocity vector u 
in two ways : 

u = Iti) + iu + ko) 

X y z 

and as 

w * J'0i + k"02 + J"03. 

Combining Equations B-10 with B-11 there results 


{UJ} .[7t]{0}, 



” 




- 


- 


- 


r 


0) 

X 

3 

cos 9 3 


-sln9 3 


sin02 





or 

0) 

y 

= 


1 



COS02 


1 


02 


0) 


sin0 3 


cos 03 



1 



L. - 


or [tt] - [T 3 ]’'[A]^. 

Now, the inverse tranformatlon of Equation B-12 is required for 
hinge kinematics applications, or it is necessary to express 
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[B-14] 


[B-15] 


[B-16] 


[B-17] 


= 


ip-1 

[A]" [T 3 ] 

-1 

([E] [E][A]'^) [T 3 ] 


= ([E][A]^) ^IE][T 3 ] 


with [E] an elementary row interchange transformation, which for 
the ( 2 - 3 - 2 ) example is 




1 


and causes 


([E][A]^) 


to be of the form; 


[E][A]^ 



such that 


([E][A]^) ^ = 

'l/a 

1 

- 


-B/a 


1 


with a » ain 02 » 
and 0 * cos 02 • 

The form of Equation B-17 is the same for all 12 types of Euler 
rotations, which was the purpose of introducing [E], and this is 
convenient with respect ot programming considerations. It follows 
that 


for types 1 , 5, 9 
for types 2 , 6 , 10 
for types 3, 7, 11 
and for types 4, 8 , 12 


a * COS 02 , 0 = sin 02 » 
a « sin02» 0 = COS02; 
a ■ -sln 02 * ^ “ COS 02 , 
a = COS 02 , 0 “ - 8 in 02 * 


Also, for each of the 12 types, there is an elementary row inter- 
change transformation [E] that can be constructed from simple 
inspection of the permutation integers of Table B-1 (2-3-2 for 
example). In fact, it is unnecessary to actually construct [E] 
because information to construct it is merely applied to [T 3 ] 
(interchanging its rows), which produces [E][T 3 ]. Thus, the 
velocity transformation of Equation B-14 can be created for any 
one of the 12 possible types with comparative ease. 
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APPENDIX C--TIME DERIVATIVES OF KINEMATIC COEFFICIENTS 


The formulation and numerical Implementation of motion equations 
for the system of interconnected bodies Involves a vector of La- 
grange multipliers, {X} (Refer to Equations II-l and II-6). In 
order to numerically evaluate {X} there is seen to be the require' 
ment of calculating time derivatives of kinematic coefficients 
(velocity transformations) associated with hinges. 

With reference to Chapter II, Section C, it is noted that for 
each hinge there is a a ^^q^ matrix of kinematic coef- 

ficients. The basic form of these matrices is repeated here, 
then the sequence of steps necessary to develop their time deriv- 
atives is Indicated. 


[C-3] 


[C-5] 


The [bp] array is 


[C-1] [bp] 


■[Tr.]“^ 

[q\] 

1 [0] 
I 

1 [1T]“^ 
1 

\q\ 

iPpj' 

[“1 

[p mj [mpj 

iw 

1 [p"”] 

i»pi 



and 


[C-2] [bq] 


"[7T]"^ 

[q^n] 

[0] 


[a. 


jp\] 


[p^n] 




Now to develop [bp] and [b^] it is necessary to expand the fol- 
lowing as : 


[A]. 

([A] [A] [p^])- 

3? ([p"»] K’]) ■ [p“»] fe'] " W [^"p']’ 


and 
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[=nf]) ■ [p"n] * [p\] ' 

The 3x3 matrix time derivatives defined by Equations C-3 through 
C“6 have factors (also 3x3 matrix time derivatives) that are ex- 
panded in terms of previously defined quantities as follows: 

[C-7] [ 6 ] - fsK* (TrI [o] {L})irRj» 

pm P m/pm 

[C-3] f R 1 • fsK* / r R 1 [a ] r 

'■ ^ qn L (qn'-q-'^njqn, 

[A] ■ ["p/0 [p'-.]’ 

with 

IC-IO) ■ SK. ([/„] ((«„) + tPpl (5„1) 

- W [/n] IV <'„>))• 

[A] ■ [pV M [p^] 

[C-12I - [si* (I\I !«„')]> 

[<=-“1 -[si* (1»<,HV)]- 

Finally, the time derivative of [r] ^ requires additional consid- 
eration. Refer to Appendix B. 

The rotation transformation [tt] is developed as 

[c-14] [7 T]“^ - |[e] [a]"^) ^ [e] [T 3 ], 

and it is shown that the form 
[C-15] ([e] [a]"^) ^ - [a] - p/a 

1 

-S/o 1 

holds for each of the 12 possible Euler rotations. In that [E] 
is constant, [A] depends only on 62 and [T 3 ] depends only on 63 , 
it follows that 
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[A] [E] [T 3 ] 

902 


[C-16] ^ [tt] ^ ■ 02 

dt 

^ k ‘^3. . 

• • 

where the Euler angle rates (62 and 63) are numerically evaluated 
before their use in Equation C-16 through application of Equation 
II-3; that is, they reside in that part of the state vector time 

derivative {y} that has been evaluated. 
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APPENDIX D—MONITOR OF SYSTEM MOMENTA AND ENERGIES 


Development of state equations for predicting dynamic response of 
a system of interconnected flexible bodies involves a consider- 
able amount of complicated formulation and programming code. This 
is certainly a true statement | independent of the particular 
method of analytical mechanics on which one might select to base 
development. 

The inherent complexity of such a digital simulation program gives 
rise to the question: is there any way of checking the program 
validity? In an attempt to answer this question, one might sug- 
gest comparing results with those of other dynamic simulations 
or hardware tests. If such a comparison is positive, then cred- 
ibility (to a degree) is established. However, there is another 
absolutely necessary (if not sufficient) condition that must be 
passed to establish validity. For a dynamic system free of ex- 
ternal forces and torques, angular and linear momenta must be 
conserved; also, total energy (kinetic plus potential) must not 
increase in time. 


It is a desirable feature for such a digital simulation program 
to have a built-in monitor of momenta and energy. The purpose of 
this appendix is to develop (in terms of previously identified 
state variables and system parameters) the expressions for total 
system angular and linear momentum vectors and the total system 
energy. 


The total angular momentum about the inertial reference can 

. 1 V 


[D-11 


expressed 
_ NB 

s-E 

j-i 


(from definition) as 


J (x X v) 


dm 


with the summation over the number of bodies (NB) of the system, 
with X being the vector positioning the elemental mass (dm) from 
the inertial origin, with v being the absolute velocity of dm, 
and with integration taken over the volume of the 

Also, from definition, the total linear momentum with 
the inertial frame is 


body 

respect to 


D-1 



[D-2] L 


NB 


V dm. 


J-1 V 


j 


Now, consistent with the notion of a body fixed axis system and 
with a consistent velocity field assumed (Refer to Chapter II, 

Section iJ.), it follows that, over the volume of the j body. 


[D-3] 


[D-4] 


[D-5] 



+ po + n , 


and 

V - + Wj x(^o +”) + ^Ic- 

on substituting Equations D—3 and D~4 Into D—1 and D— 2 and Integra' 
ting, it becomes clear that the first six elements of the product 


{p}j =■ Mj 


{Pv> 

{P|} 


are projections of the body's angular and linear momentum 
vectors onto the moving body axis system. In fact, tP^^J includes 

the effect of momentum wheels (See Equation 11—109), which surely 
must be accounted for. 


Thus, the angular momentum of the body (about its body— origin) 

is 


[D-6] -nj - [.J 

while the linear momentum of the body is 


[D-71 

where J is the unit vector basis associated with the body fixed 
reference triad. 


D-2 



¥ 


Now rotation transformations that relate vector components in 
each body system to the inertial system exist; also, position 
vector from the inertial origin to the reference point of each 
body exists. It follows that 
_ NB _ 

NB 

’1} M ‘'■''’r 

and that 

_ NB . 

[D-9] 

J-1 

NB , 

■ £ ( [°“j] ‘““’j " ["''* ( ["“J Wj) • 


The total angular and linear momentum vectors are calculated by 
the program in the manner indicated in Equations D-8 and D-9. For 
a variety of torque/force-free configurations that have been ex- 
amined, momentum has been conserved within acceptable numerical 
tolerances. 


The total energy is calculated (Refer to Equations 11-38 and II- 
42.) as 

[D-10] T + V = I [yjj [m]j jujj + 

j=l 

[sjj Wj (Uj) . 

The kinetic energy contribution of embedded momentum wheels is 
included (as it must be), because Includes momentum wheel 

inertial coupling terms and Includes momentum wheel spin 

rates (5 ). 
s 

Potential energy, additional to that shown in Equation D-10, comes 
about in the event that there is a **sprung’* hinge; say for example, 
associated with the coordinate. If the spring force/torque is 

linear with 6^^, then additional potential energy is 

[D-ll] \/ - K, . 

^ (additional) k k I >.3 



APPENDIX E—ROOT LOCUS SOLUTION TECHNIQUES 


The Root Locus solution procedure that is included in the digital 
program was supplied by Goddard Space Flight Center (GSFC) per- 
sonnel and is included here for completeness. The following dis- 
cussion has been extracted from GSFC Branch Report No 254 dated 
October 2, 1973 (Mr. Raymond V. Welch, author). 


Two expressions are identified as the basis for initiating the 
solution process. 


[E-l] 


P(s) 

Q(s) 


1 


[E-2] 



180° K > 0 


Expressions E-l and E-2 contain polynominal expressions for the 
conventional open loop expression 


[E-3] KG(s) H(s) » KP(s)/Q(s). 


The solution process requires a starting point (known to lie 
somewhere on the loci) from which to generate the desired locus; 
therefore, assume there exists a known value of s, say s , such 
that ° 


[E-4] /pCSq) 

i.e., s^ is a point on the locus. A good starting point might be 

an open loop pole or zero as these points are usually known a 
priori; however, any point along the branch of the locus to be 
determined may be used. The locus is then traced using the fol- 
lowing procedure. 


Step 1 

Draw a "small** circle of radius r centered at s and define s 

o c 

to be the values of s which lie on this circle (See Figure E-l.). 
These values of s are determined analytically by: 

[E-5] s = s + 

^ ^ c 0 


E-l 


where 6 Is meausred as the angle generated by a counterclockwise 
rotation from the positive real axis of the s^plane. If the locus 
does not terminate inside the circle, then there must exist at 
least two values of 0 between zero and 360 degrees such that 


lE-6] 



Q(Sc) 


180® 


(More than two values of 0 will exist that satisfy this equation 
if breakway points of the locus are encircled or if the circle is 
large enough to Intersect other branches of the locus. Consequently, 
to avoid changing branches, r must be kept smaller than the distance 
between the branches of the locus). Suppose that for 0 - and 

0 * 02 the above equation is satisfied, then s _ and s « are points 

cl c2 

on the locus where 


[E“7] s - = 3 + re^^^ 

cl o 

s - = s + re ^ 
c2 o 

See Figure E-1. These roots are found by an iterative process. 
Step 2 

Define \fj(6) by 


[E-8] 


^( 0 ) 




where s^ is defined above and 0 is any arbitrary angle. If ^(0) 

does not equal 180 degrees increment 0 by A9 and reevaluate 
Continue this process until ^(0) passes through 180 degrees. 

(ij;(0) is a monotonic function across the locus, thus ij/(0) crossing 
180 degrees implies the locus has been crossed). When ip(0) passes 
through 180 degrees, redefine A0 as 


[E-9] A0 » -KA0 - 0<K<1 

^ new old 

and again calculate iK9). If ^(0) does not equal 180 degrees for 
this point, increment 0 by this new A9 and recalculate ^(9). Con- 
tinue this process until ip(0) again crosses 180 degrees. Reduce A9 
again as above and repeat the previous operation until A0 •= 6 , 
where 5 is some predefined small positive number. At this time 
t/;( 0) will be 180 ±e degrees, where c is a "small” angle whose size 
is a function of 6; thus a root on the locus has been found. This 
root is either s^^ or depending on the initial choice of 0 

and on the initial sign of A0 . For this subprogram, the initial 



Figure E~1 Variable Definitions at Starting Point of Locus 
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Figure E-2 Search Area Definition for Finding Roots After the 
First Root is Found. 

E-4 


r* 


[E-10] 


[E-11] 


[E-12] 


value of 0 la optional, but if no choice is made a value of 180 
degrees is used. Also the initial increment A0 is set at 10 de- 
grees with the change in A0 for every 180 degrees crossing defined 
by 


A0 

new 



old 


It was found that 6 = 10^® degrees was sufficiently small to 
Insure e< .001 degrees except possibly for values of s near the 
open loop poles or zeros. Assume that 0 is chosen initially so 
that the root found is located at the point s - as defined above 
. and as shown in Figure E-1, and suppose it is desired to continue 
the locus in this direction. 


Step 3 


Draw a circle of radius _r centered at Again if the locus 

docs not terminate inside the circle, the locus will intersect 

the circle in at least two points. Because the circles located 

at s and s . are of equal radius, one of these points is s as 
o cl o 

shown in Figure E-2. This root is eliminated from the search 

routine by restricting the range of 0 to 

01 -120° <0 <01 -fl20° 

These limits are chosen as they are the points where the circles 

located at s and s , intersect. Thus the search for a new root 
o cl 

is conducted only in a previously unsearched region. The initial 
angle is chosen at 0 ■ 0i - 120 degrees with A6>o and step 2 is 
repeated until A0 » 6, i.e., another root is found. 


Step 4 

Draw a circle of radius r centered at the root found in the pre- 
vious step and repeat Step 2 with the restrictions on 0 as defined 
in Step 3. 


Repeating Step 4, the roots along one branch of the locus are de- 
termined. The value of the gain, K, for each of these roots is 
calculated from 


K 


Q(s) 
P(s) • 


E-5 
















